{ "cells": [ { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "view-in-github" }, "source": [ "\"Open" ] }, { "cell_type": "markdown", "metadata": { "id": "tT1kZRa2aNbd" }, "source": [ "# 2 Step Adam Bashforth \n", "#### John S Butler \n", "john.s.butler@tudublin.ie \n", "[Course Notes](https://johnsbutler.netlify.com/files/Teaching/Numerical_Analysis_for_Differential_Equations.pdf) [Github](https://github.com/john-s-butler-dit/Numerical-Analysis-Python)\n", "\n", "\n", "\n", "This notebook implements the 2 step Adams Bashforth method for three different population intial value problems.\n", "\n", "# Formula\n", "The general 2 step Adams-Bashforth method for the first order differential equation\n", "\\begin{equation} y^{'} = f(t,y), \\end{equation}\n", "numerical approximates $y$ the at time point $t_i$ as $w_i$\n", "with the formula:\n", "\\begin{equation} w_{i+1}=w_i+\\frac{h}{2}\\big[3f(t_i,w_i)-f(t_{i-1},w_{i-1})\\big],\\end{equation}\n", "for $i=0,...,N-1$, where \n", "\n", "and $h$ is the stepsize.\n", "\n", "To illustrate the method we will apply it to three intial value problems:\n", "## 1. Linear \n", "Consider the linear population Differential Equation\n", "\\begin{equation} y^{'}=0.1y, \\ \\ (2000 \\leq t \\leq 2020), \\end{equation}\n", "with the initial condition,\n", "\\begin{equation}y(2000)=6.\\end{equation}\n", "\n", "## 2. Non-Linear Population Equation \n", "Consider the non-linear population Differential Equation\n", "\\begin{equation} y^{'}=0.2y-0.01y^2, \\ \\ (2000 \\leq t \\leq 2020), \\end{equation}\n", "with the initial condition,\n", "\\begin{equation}y(2000)=6.\\end{equation}\n", "\n", "## 3. Non-Linear Population Equation with an oscillation \n", "Consider the non-linear population Differential Equation with an oscillation \n", "\\begin{equation} y^{'}=0.2y-0.01y^2+\\sin(2\\pi t), \\ \\ (2000 \\leq t \\leq 2020), \\end{equation}\n", "with the initial condition,\n", "\\begin{equation}y(2000)=6.\\end{equation}" ] }, { "cell_type": "markdown", "metadata": { "id": "Ea_Jx1iBaNbf" }, "source": [ "#### Setting up Libraries" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "id": "cL3QkPThaNbg" }, "outputs": [], "source": [ "## Library\n", "import numpy as np\n", "import math \n", "import pandas as pd\n", "\n", "%matplotlib inline\n", "import matplotlib.pyplot as plt # side-stepping mpl backend\n", "import matplotlib.gridspec as gridspec # subplots\n", "import warnings\n", "\n", "warnings.filterwarnings(\"ignore\")" ] }, { "cell_type": "markdown", "metadata": { "id": "SYUkyLUiaNbl" }, "source": [ "## Discrete Interval\n", "The continuous time $a\\leq t \\leq b $ is discretised into $N$ points seperated by a constant stepsize\n", "\\begin{equation} h=\\frac{b-a}{N}.\\end{equation}\n", "Here the interval is $2000\\leq t \\leq 2020,$ \n", "\\begin{equation} h=\\frac{2020-2000}{200}=0.1.\\end{equation}\n", "This gives the 201 discrete points:\n", "\\begin{equation} t_0=2000, \\ t_1=2000.1, \\ ... t_{200}=2020. \\end{equation}\n", "This is generalised to \n", "\\begin{equation} t_i=2000+i0.1, \\ \\ \\ i=0,1,...,200.\\end{equation}\n", "The plot below shows the discrete time steps:" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 298 }, "id": "gf-GKCD4aNbm", "outputId": "e5ce652a-6de4-47bb-fecb-8c9ac5768dc6" }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmEAAAEICAYAAAAX5iNEAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAfh0lEQVR4nO3deZxkZX3v8c8XUBBRGBYRBGZcSAzGvV0w3ohREUwQY0iCMTqokXCJuZrEJBhjNLhrosRIosQlBLmC0RgnMVzEheReVLRREFGRJcriACOrgwqiv/vHOS01TVUv0zX9dPd83q9XveYsTz1LPae6vn3OqZ5UFZIkSVpc27TugCRJ0tbIECZJktSAIUySJKkBQ5gkSVIDhjBJkqQGDGGSJEkNGMKkXpKjkvy/gfVK8qCWfRolybuSvKpBu/8zybVJNibZbQ7lv5Xkqf3ynyV5z5bv5ZaT5LlJPtG6HzNJsl8/P9suQlu/muTKvr1HjqG+s5P8zjj6Ji0HhjBt1QZDwhaq/x+TvG6BdWwSDgGq6piqeu3CejfvftwNeBtwcFXtVFXXz+f5VfWGqlq0D9gka/ogvd24nl9Vp1bVwePr5fhV1RX9/Px4trILfY2AvwJe0rf35c2sY+ySnJTk4iQ/SXLULGW3T/K+JLckuSbJHy5SNyVDmNTSAj78WtgT2AG4qHVHYNm9divVajbzeNjCZ+ouAI4FvjSHsq8B9qcby5OBP0lyyJbrmnQnQ5g0B9MvkwyenUrn7Umu63+bvjDJzyc5Gngu3Q/1jUn+rS//rSR/muQrwK1JtktyXJLLknwvydeS/Gpf9ueAdwEH9nXc1G/f5AxbkhcnuTTJDUnWJdl7YF8lOSbJJUluSnJikowY5/ZJTkjynf5xQr/tZ4CL+2I3Jfn0iOc/L8m3k1yf5JXT9r0myQf65R2SfKAvd1OSLybZs9+3a5L39+3fmORf++0HJbmqf+2uAd6fZJuB1+76JB9Ksmvf5H8N9HdjkgP7el6Y5Ot93WcmWT1i2u/y/OlnJfvX9tj+tf1ektcmeWCSz/bHwoeS3H2g/K8kOb8f82eTPGxE21N1/68klyf5bpK3Jtmm37dNkj/vX+vrkvxTkp37fZuc3eqP3dcmOafv4yeS7D7DGB+U5D+T3Ny3e/qQvm2fZCOwLXBBksv67T/Xt3dTkouSPHPgOf+Y5O+T/EeSW+kCzzCrR/R1zqrqxKr6FPDDORRfC7y2qm6sqq8D/wAcNd82pc1hCJMW7mDgF4GfAXYGfgO4vqpOAk4F3tJfrjls4DnPAX4Z2KWq7gAuA/5H//y/BD6QZK/+Q+EY4HN9HbtMbzzJLwFv7NvdC/g2cNq0Yr8CPAZ4WF/u6SPG8krg8cAjgIcDjwX+vKq+CTykL7NLVf3SkH4cAPw98Dxgb2A3YJ8R7aztx7pvX+4Y4Af9vlOAHfv27gO8feB59wV2pTtrcTTw+8CzgCf1bd4InNiX/cWB/u5UVZ9LcjjwZ8CzgT2A/wt8cEQf7/L8EeWeDjya7nX7E+Ak4Lf7sf083VyT7p6p9wG/24/53cC6JNuPqBfgV4EJ4FHA4cAL++1H9Y8nAw8AdgLeOUM9vwW8gO71vDvw8hnG+FrgE8Aquvn72+mVVdVtVbVTv/rwqnpgusvV/9Y/9z50c3Nqkp+d1o/XA/cCNrnEPoe+0oe7UY/jZhj/UElW0b1nLhjYfAF3HuvSFmUIkxbuR3QfKg8GUlVfr6r1szznHVV1ZVX9AKCq/rmqvlNVP6mq04FL6ALQXDwXeF9VfamqbgNeQXfmbM1AmTdV1U1VdQXwGbqQNaqu46vquqraQBcInzfHfhwB/HtV/Vffj1cBPxlR9kd0QeRBVfXjqjqvqm5JshdwKHBMf2biR1X1nwPP+wnw6j4E/IAuvL2yqq7q23wNcERGX6o8BnhjP0d3AG8AHjHD2bC5eEtV3VJVFwFfBT5RVZdX1c3AGcDUDetHA++uqnP7MZ8M3EYX3kZ5c1Xd0M/bCfSBjm6e3ta3s5Fuzo+cYdzvr6pv9q/Zhxg9/9DNzWpg76r6YVWNCkvTPZ4uDL6pqm6vqk8D/z7QZ4CPVdU5/XE+6izVyL5W1S4zPN40x34OmgqSNw9su5nu/SxtcYYwaYH6D5t30p2BuS7dTcH3nuVpVw6uJHn+wGWqm+jOoMz1MszedGe/pvqzEbgeuN9AmWsGlr/PnR8+M9bVL+89ouyw5/50XFV1a9+PYU4BzgRO6y87vqU/k7IvcENV3TjieRumfXivBj468Lp9Hfgx3f1rw6wG/mag/A1A2PS1mq9rB5Z/MGR96rVeDfzR4NkbuvHO9PoOHieDczFsnrZj9LjnOv/Qnc0L8IX+kuILZyg7aG/gyqoaDN7fZtPX9kpmN5++LtTG/t/B9+u9ge9twTalnzKESXNzK90lsin3HdxZVe+oqkcDB9BdlvzjqV0j6vvp9v4szD8ALwF26y85fpXug3CmOqZ8h+4Dfqq+e9KdZbp6lufNWhewX79tLtbThYqpfuzY9+Mu+jNcf1lVBwBPoLtc+ny6D+ldk9zlsuvUU6etXwkcOu2MyA5VdfWQslPlf3da+XtU1Wfn0NZCXQm8flrbO1bVqMuhMPB6sulcDJunO9g0AM7FXcZYVddU1Yuram+6S6d/l7n9qZbvAPtO3bc20K/B43BBr2l/39qox5/Nt74+7K+nu/Q+5eEskS+faOUzhElzcz7w7CQ79h9IL5rakeQxSR7Xn8m5le5m4KmzAdfS3bMzk3vSfTht6Ot7Ad2ZsCnXAvsM3uA9zQeBFyR5RH9/0RuAc6vqW/MZ4EBdf55kj/6G6L8APjDH534Y+JUkT+z7ejwjfsYkeXKSh6b7htwtdJfAftJfxj2D7oN/VZK7JfnFYXX03gW8fupyYt/vw/t9G+jm4QHTyr8iyUP68jsn+fURdQ97/kL8A3BMf6wkyT2T/HKSmS59/XH/OuwLvBSYukn+g8AfJLl/kp3o5vz0/hLrfNxljEl+PcnUvXw30h2boy4rDzqX7szVn/TzdhBwGHe9P3Gz9fetjXq8YWAMd0+yA90vMndL90WQUZ93/0R3zK9K8mDgxcA/jqvP0kwMYdLcvB24nS4QnUx3w/2Ue9N9wN5Id/nleuCt/b73Agf0l5/+dVjFVfU14K+Bz/X1PxQ4Z6DIp+l+M78myXeHPP+TdPdffYTut/oHAkdu1ijhdcAk8BXgQrqv+M/p75z190T9HvC/+37cCFw1ovh96ULbLXSXEP+T7hIldPeg/Qj4BnAd8LIZmv0bYB3wiSTfAz4PPK7vz/fpbgI/p3/9H19VHwXeTHcZ9Ba6M46HjhjPXZ4/64swg6qapPuAfyfda3Mps38L72PAeXS/BHyc7niC7gb/U+i+3fjfdMH/9zejT8PG+Bjg3HTfflwHvLSqLp9DXbfTha5Dge8Cfwc8v6q+Md9+jcEn6C4FP4HuixI/oP8SQro/uDt4puvVdF+M+TbdcfjWqvo/i9tdba1SNe4z7pKkhUpSwP5VdWnrvkjaMjwTJkmS1IAhTJIkqQEvR0qSJDXgmTBJkqQGluV/gLv77rvXmjVrWndDkiRpVuedd953q2qP6duXZQhbs2YNk5OTrbshSZI0qyTfHrbdy5GSJEkNGMIkSZIaMIRJkiQ1YAiTJElqwBAmSZLUgCFMkiSpAUOYJElSA4YwSZKkBgxhkiRJDRjCJEmSGjCESZIkNWAIkyRJasAQJkmS1IAhTJIkqQFDmCRJUgOGMEmSpAYMYZIkSQ0YwiRJkhowhEmSJDVgCJMkSWrAECZJktSAIUySJKkBQ5gkSVIDhjBJkqQGDGGSJEkNjCWEJTkkycVJLk1y3JD92yc5vd9/bpI10/bvl2RjkpePoz+SJElL3YJDWJJtgROBQ4EDgOckOWBasRcBN1bVg4C3A2+etv9twBkL7YskSdJyMY4zYY8FLq2qy6vqduA04PBpZQ4HTu6XPww8JUkAkjwL+G/gojH0RZIkaVkYRwi7H3DlwPpV/bahZarqDuBmYLckOwF/CvzlbI0kOTrJZJLJDRs2jKHbkiRJ7bS+Mf81wNurauNsBavqpKqaqKqJPfbYY8v3TJIkaQvabgx1XA3sO7C+T79tWJmrkmwH7AxcDzwOOCLJW4BdgJ8k+WFVvXMM/ZIkSVqyxhHCvgjsn+T+dGHrSOC3ppVZB6wFPgccAXy6qgr4H1MFkrwG2GgAkyRJW4MFh7CquiPJS4AzgW2B91XVRUmOByarah3wXuCUJJcCN9AFNUmSpK1WuhNSy8vExERNTk627oYkSdKskpxXVRPTt7e+MV+SJGmrZAiTJElqwBAmSZLUgCFMkiSpAUOYJElSA4YwSZKkBgxhkiRJDRjCJEmSGjCESZIkNWAIkyRJasAQJkmS1IAhTJIkqQFDmCRJUgOGMEmSpAYMYZIkSQ0YwiRJkhowhEmSJDVgCJMkSWrAECZJktSAIUySJKkBQ5gkSVIDhjBJkqQGDGGSJEkNGMIkSZIaMIRJkiQ1YAiTJElqwBAmSZLUgCFMkiSpAUOYJElSA4YwSZKkBsYSwpIckuTiJJcmOW7I/u2TnN7vPzfJmn7705Kcl+TC/t9fGkd/JEmSlroFh7Ak2wInAocCBwDPSXLAtGIvAm6sqgcBbwfe3G//LnBYVT0UWAucstD+SJIkLQfjOBP2WODSqrq8qm4HTgMOn1bmcODkfvnDwFOSpKq+XFXf6bdfBNwjyfZj6JMkSdKSNo4Qdj/gyoH1q/ptQ8tU1R3AzcBu08r8GvClqrptDH2SJEla0rZr3QGAJA+hu0R58AxljgaOBthvv/0WqWeSJElbxjjOhF0N7Duwvk+/bWiZJNsBOwPX9+v7AB8Fnl9Vl41qpKpOqqqJqprYY489xtBtSZKkdsYRwr4I7J/k/knuDhwJrJtWZh3djfcARwCfrqpKsgvwceC4qjpnDH2RJElaFhYcwvp7vF4CnAl8HfhQVV2U5Pgkz+yLvRfYLcmlwB8CU3/G4iXAg4C/SHJ+/7jPQvskSZK01KWqWvdh3iYmJmpycrJ1NyRJkmaV5Lyqmpi+3b+YL0mS1IAhTJIkqQFDmCRJUgOGMEmSpAYMYZIkSQ0YwiRJkhowhEmSJDVgCJMkSWrAECZJktSAIUySJKkBQ5gkSVIDhjBJkqQGDGGSJEkNGMIkSZIaMIRJkiQ1YAiTJElqwBAmSZLUgCFMkiSpAUOYJElSA4YwSZKkBgxhkiRJDRjCJEmSGjCESZIkNWAIkyRJasAQJkmS1IAhTJIkqQFDmCRJUgOGMEmSpAYMYZIkSQ0YwiRJkhowhEmSJDUwlhCW5JAkFye5NMlxQ/Zvn+T0fv+5SdYM7HtFv/3iJE8fR38W5NRTYc0a2Gab7t9TT12ebSxWO45l6bWxWO2slDYWq52V0sZiteNYll4bi9XOSmljLqpqQQ9gW+Ay4AHA3YELgAOmlTkWeFe/fCRwer98QF9+e+D+fT3bztbmox/96NoiPvCBqh13rII7Hzvu2G1fTm0sVjuOZem1sVjtrJQ2FqudldLGYrXjWJZeG4vVzkppYxpgsoZlqGEb5/MADgTOHFh/BfCKaWXOBA7sl7cDvgtketnBcjM9tlgIW71600mZeuy3X9WTnlR1yilduVtv7dZPO61bv+mmbv0jH+nWN2zo1tet69bXr+/WzzhjdBvbb1919tld+W98oyt/zjnd+oUXdutf+EK3/uUvd+tf/nK3/oUvdOsXXtitn3NOV9+wdlavrjrrrK78FVd05c84o1tfv75bX7euW9+woVv/yEe69Ztu6tZPO61b33ff4W3stlu3f8pJJ1U95Sl3rp94YtUhh9y5fsIJVYcdduf6W99a9exn37m+yy6jx1JV9apXVR111J3ljzuu6sUvvnP9j/6o6thj71x/6Uu7x5Rjj626171Gt3HUUV0bU5773Krjj79z/Td/s+qNb7xz/dnP7sYw5bDDujGOmvsdduheoylPelLV+9/fLd9++/yPvZnm/ooruvJnndWVv+yybn2+x95eew1vY889u/3f+EZX/uyzu/XLLuvW53PszfReufXWrvwpp3Tlb7+9W3//++d/7N3jHqNfrze+sZvfKccf383/lLkee6PGstNOXZkpL35xV8eU+Rx7o9pYterO8occ0r0GU57ylPkfe7vvPrydffYZ/XOvan7H3kxzP9PPvfkeezO9V0b93JvvsTfTvMz0c28+x95Mx9eUYT/35nvs7bzz6Ndr1M+9KXM99maa+835zK2667G3994zf6ZsAaNC2DguR94PuHJg/ap+29AyVXUHcDOw2xyfC0CSo5NMJpncsGHDGLo9xBVXDN9+5ZXDt4+zjdtuG18bM9U3qv3NcdVVw7dff/342gC46abh28c5lu99b8u3MaquH/5wfG3A4sz9NdcM337ddeNrY7HeKz/4wfza3xyj6tq4ccu3ceON42sDRr+/r756fG2spJ+TizEvi3F8Adx88/za3xyLMffr18+v7S1pWDKbzwM4AnjPwPrzgHdOK/NVYJ+B9cuA3YF3Ar89sP29wBGztbnoZ8LGmY4Xo43FasexLL02FqudldLGYrWzUtpYrHYcy9JrY7HaWSltTMMWPBN2NbDvwPo+/bahZZJsB+wMXD/H5y6e178edtxx02077thtX05tLFY7jmXptbFY7ayUNharnZXSxmK141iWXhuL1c5KaWOuhiWz+Tzo7vG6nO7G+qkb8x8yrczvsemN+R/qlx/CpjfmX07LG/OruhvzVq+uSrp/t8SNeovRxmK141iWXhuL1c5KaWOx2lkpbSxWO45l6bWxWO2slDYGMOJMWLp9C5PkGcAJdN+UfF9VvT7J8X2j65LsAJwCPBK4ATiyqi7vn/tK4IXAHcDLquqM2dqbmJioycnJBfdbkiRpS0tyXlVN3GX7OELYYjOESZKk5WJUCPMv5kuSJDVgCJMkSWrAECZJktSAIUySJKkBQ5gkSVIDhjBJkqQGDGGSJEkNGMIkSZIaMIRJkiQ1YAiTJElqwBAmSZLUgCFMkiSpAUOYJElSA4YwSZKkBgxhkiRJDRjCJEmSGjCESZIkNWAIkyRJasAQJkmS1IAhTJIkqQFDmCRJUgOGMEmSpAYMYZIkSQ0YwiRJkhowhEmSJDVgCJMkSWrAECZJktSAIUySJKkBQ5gkSVIDhjBJkqQGFhTCkuya5Kwkl/T/rhpRbm1f5pIka/ttOyb5eJJvJLkoyZsW0hdJkqTlZKFnwo4DPlVV+wOf6tc3kWRX4NXA44DHAq8eCGt/VVUPBh4J/EKSQxfYH0mSpGVhoSHscODkfvlk4FlDyjwdOKuqbqiqG4GzgEOq6vtV9RmAqrod+BKwzwL7I0mStCwsNITtWVXr++VrgD2HlLkfcOXA+lX9tp9KsgtwGN3ZNEmSpBVvu9kKJPkkcN8hu145uFJVlaTm24Ek2wEfBN5RVZfPUO5o4GiA/fbbb77NSJIkLSmzhrCqeuqofUmuTbJXVa1Pshdw3ZBiVwMHDazvA5w9sH4ScElVnTBLP07qyzIxMTHvsCdJkrSULPRy5Dpgbb+8FvjYkDJnAgcnWdXfkH9wv40krwN2Bl62wH5IkiQtKwsNYW8CnpbkEuCp/TpJJpK8B6CqbgBeC3yxfxxfVTck2YfukuYBwJeSnJ/kdxbYH0mSpGUhVcvvyt7ExERNTk627oYkSdKskpxXVRPTt/sX8yVJkhowhEmSJDVgCJMkSWrAECZJktSAIUySJKkBQ5gkSVIDhjBJkqQGDGGSJEkNGMIkSZIaMIRJkiQ1YAiTJElqwBAmSZLUgCFMkiSpAUOYJElSA4YwSZKkBgxhkiRJDRjCJEmSGjCESZIkNWAIkyRJasAQJkmS1IAhTJIkqQFDmCRJUgOGMEmSpAYMYZIkSQ0YwiRJkhowhEmSJDVgCJMkSWrAECZJktSAIUySJKkBQ5gkSVIDhjBJkqQGFhTCkuya5Kwkl/T/rhpRbm1f5pIka4fsX5fkqwvpiyRJ0nKy0DNhxwGfqqr9gU/165tIsivwauBxwGOBVw+GtSTPBjYusB+SJEnLykJD2OHAyf3yycCzhpR5OnBWVd1QVTcCZwGHACTZCfhD4HUL7IckSdKystAQtmdVre+XrwH2HFLmfsCVA+tX9dsAXgv8NfD92RpKcnSSySSTGzZsWECXJUmS2ttutgJJPgncd8iuVw6uVFUlqbk2nOQRwAOr6g+SrJmtfFWdBJwEMDExMed2JEmSlqJZQ1hVPXXUviTXJtmrqtYn2Qu4bkixq4GDBtb3Ac4GDgQmknyr78d9kpxdVQchSZK0wi30cuQ6YOrbjmuBjw0pcyZwcJJV/Q35BwNnVtXfV9XeVbUGeCLwTQOYJEnaWiw0hL0JeFqSS4Cn9uskmUjyHoCquoHu3q8v9o/j+22SJElbrVQtv9urJiYmanJysnU3JEmSZpXkvKqamL7dv5gvSZLUgCFMkiSpAUOYJElSA4YwSZKkBgxhkiRJDRjCJEmSGjCESZIkNWAIkyRJasAQJkmS1IAhTJIkqQFDmCRJUgOGMEmSpAYMYZIkSQ0YwiRJkhowhEmSJDVgCJMkSWrAECZJktSAIUySJKkBQ5gkSVIDhjBJkqQGDGGSJEkNGMIkSZIaMIRJkiQ1YAiTJElqIFXVug/zlmQD8O0t3MzuwHe3cBtL1dY8dti6x781jx227vFvzWOHrXv8jn3LW11Ve0zfuCxD2GJIMllVE6370cLWPHbYuse/NY8dtu7xb81jh617/I693di9HClJktSAIUySJKkBQ9hoJ7XuQENb89hh6x7/1jx22LrHvzWPHbbu8Tv2RrwnTJIkqQHPhEmSJDVgCJMkSWpgxYawJPsm+UySryW5KMlL++27JjkrySX9v6v67UnyjiSXJvlKkkcN1LW2L39JkrUj2htabwvjGnuSRyT5XF/HV5L85oj2jkqyIcn5/eN3Fm+0d+nLOOf9xwNjWjeive2TnN4//9wkaxZjnKOMce6fPDD285P8MMmzhrS3nOf+wf3xfVuSl0+r65AkF/evy3Ej2lsycz+usY+qZ0h7ByW5eWDe/2JxRjrcmOf+W0ku7Mc1OaK9kT83FtsY5/5np73nb0nysiHtLfe5f24/Zxcm+WyShw/Utfjv+6pakQ9gL+BR/fK9gG8CBwBvAY7rtx8HvLlffgZwBhDg8cC5/fZdgcv7f1f1y6uGtDe03mU+9p8B9u+X9wbWA7sMae8o4J2t53ycY+/3bZxDe8cC7+qXjwROXynjH6hzV+AGYMcVNvf3AR4DvB54+UA92wKXAQ8A7g5cABywlOd+jGMfWs+Q9g4C/r31nI97/P2+bwG7z9LerO+b5Tj2gTq3Ba6h+wOjK23un0D/GQ4cyp2fd03e981fwEWcqI8BTwMuBvYamLyL++V3A88ZKH9xv/85wLsHtm9Sbnr56fUuhcfmjn1IPRfQh7Jp249iiXwQj3PszC2EnQkc2C9vR/eXl9N63OOce+Bo4NQR9S/buR8o9xo2DSIHAmcOrL8CeMVymvvNHfuoeoZsP4gl9EE8zvEztxA2p5+Zy23sA/sOBs4ZsW9FzH2/fRVwdb/c5H2/Yi9HDupPFz4SOBfYs6rW97uuAfbsl+8HXDnwtKv6baO2Tzeq3qYWOPbBeh5L99vBZSOa+rX+FO+Hk+w7nt4vzBjGvkOSySSfz5BLcdOfX1V3ADcDu41rDAsxrrmn+23vgzM0tVznfpS5vueX5NwvcOyj6hnmwCQXJDkjyUM2t7/jNobxF/CJJOclOXpEmbkeI4tqXHPP7O/5lTL3L6I7owmN3vcrPoQl2Qn4CPCyqrplcF91UXbsf6NjS9U7X+Mae5K9gFOAF1TVT4YU+TdgTVU9DDgLOHlBHR+DMY19dXX/ncVvASckeeD4e7pljHnuH0r3298wK3Xul6UxzvvIenpfont/PBz4W+BfF9TxMRnT+J9YVY+iu1T1e0l+cfw9Hb8xzv3dgWcC/zyiyIqY+yRPpgthf7ponRxiRYewJHejm5RTq+pf+s3X9h8sUx8w1/XbrwYGf4vfp982avt0o+ptYkxjJ8m9gY8Dr6yqzw9rq6qur6rb+tX3AI8e51jma1xjr6qpfy8Hzqb7DWu6nz4/yXbAzsD1YxzOvI1r/L3fAD5aVT8a1tYyn/tR5vqeX1JzP6axj6pnE1V1S1Vt7Jf/A7hbkt3HMIzNNq7xD7zvrwM+Cjx2SLG5HiOLYlxj7x0KfKmqrh22cyXMfZKH0f28Oryqpt6zTd73KzaEJQnwXuDrVfW2gV3rgLX98lq668dT25+fzuOBm/tTmWcCBydZ1X+74mCGnxUYVe+iG9fY+9+IPgr8U1V9eIb29hpYfSbw9TENZd7GOPZVSbbv69wd+AXga0OaHKz3CODT/W9dTYzxuJ/yHGa4LLHM536ULwL7J7l//x44sq9juiUz9+Ma+wz1TC93377s1K0K29A2gI5r/PdMcq+pZbqf918dUnS2982iGeNxP2W29/yynvsk+wH/Ajyvqr45UL7N+35zbyZb6g/giXSnH78CnN8/nkF37fZTwCXAJ4Fd+/IBTqS75+lCYGKgrhcCl/aPFwxsf89UuVH1LuexA78N/GigjvOBR/T7jgee2S+/EbiI7sb9zwAPXgFjf0K/fkH/74sG2hgc+w50p+0vBb4APGAFHfdr6H7r22ZaGytl7u9Ld9/HLcBN/fK9+33PoPuW1WV0Z4GX9NyPa+yj6umfcwxwTL/8koF5/zzwhGV23I8a/wP6MV3Qj29w7gfHP/J9s1zH3u+7J12g2nlaGytp7t8D3DhQdnKgrkV/3/vfFkmSJDWwYi9HSpIkLWWGMEmSpAYMYZIkSQ0YwiRJkhowhEmSJDVgCJMkSWrAECZJktTA/wcq1l6BcGAlrAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" }, { "data": { "text/plain": [ "21" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "### DISCRETE TIME\n", "N=20\n", "t_end=2020.0\n", "t_start=2000.0\n", "h=((t_end-t_start)/N)\n", "t=np.arange(t_start,t_end+h/2,h)\n", "\n", "## PLOTS TIME\n", "fig = plt.figure(figsize=(10,4))\n", "plt.plot(t,0*t,'o:',color='red')\n", "plt.title('Illustration of discrete time points for h=%s'%(h))\n", "plt.show()\n", "len(t)" ] }, { "cell_type": "markdown", "metadata": { "id": "v28phthjaNbq" }, "source": [ "# 1. Linear Population Equation\n", "## Exact Solution \n", "The linear population equation\n", "\\begin{equation} y^{'}=0.1y, \\ \\ (2000 \\leq t \\leq 2020), \\end{equation}\n", "with the initial condition,\n", "\\begin{equation}y(2000)=6.\\end{equation}\n", "has a known exact (analytic) solution\n", "\\begin{equation} y(t)=6e^{0.1(t-2000)}. \\end{equation}\n", "\n", "## Specific 2 step Adams Bashforth\n", "The specific 2 step Adams Bashforth for the linear population equation is:\n", "\n", "\\begin{equation}w_{i+1}=w_{i}+\\frac{h}{2}\\big[3(f(t_i,w_i))-f(t_{i-1},w_{i-1}))\\big] \\end{equation}\n", "where\n", "\\begin{equation}f(t,y)=0.1y.\\end{equation}" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "id": "gY78y12QaNbq" }, "outputs": [], "source": [ "## THIS IS THE RIGHT HANDSIDE OF THE LINEAR POPULATION DIFFERENTIAL \n", "## EQUATION\n", "def linfun(t,w):\n", " ftw=0.1*w\n", " return ftw" ] }, { "cell_type": "markdown", "metadata": { "id": "nWsedHwUaNbs" }, "source": [ "this gives\n", "\n", "\\begin{equation} w_{i+1}=w_{i}+\\frac{0.1}{2}\\big[ 3(0.1w_i)-0.1w_{i-1} \\big]\\end{equation}\n", "for $i=1,...,199$, where $w_i$ is the numerical approximation of $y$ at time $t_i$, with step size $h$ and the initial condition\n", "\\begin{equation}w_0=6.\\end{equation}" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "id": "JluvNM93aNbt" }, "outputs": [], "source": [ "### INSERT METHOD HERE\n", "w=np.zeros(N+1) # a list of 2000+1 zeros\n", "w[0]=6 # INITIAL CONDITION\n", "w[1]=6.06\n", "for i in range(1,N):\n", " w[i+1]=w[i]+h/2*(3*linfun(t[i],w[i])-linfun(t[i-1],w[i-1]))\n" ] }, { "cell_type": "markdown", "metadata": { "id": "UXkHGtmuaNbv" }, "source": [ "## Plotting Results" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 265 }, "id": "kzfRp2p1aNbv", "outputId": "6bac5104-2832-4707-ee6c-c6325966e480" }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeAAAAD4CAYAAAA0JjXXAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXxU1f3/8dfJJBCysgSSSFgiyqaRRUARBRHRVgkoti5gEVz4un61rVYtv7Zu6Fdbt/r9ao0LWMRC3UBcSyOKGwgqEEgM+05CwBASIPv5/TGTNIGZZJLcJJPk/Xw85sHMnTufe05uwmfuuWcx1lpERESkaQU1dwFERETaIiVgERGRZqAELCIi0gyUgEVERJqBErCIiEgzCG7Kg8XExNjevXs35SFFRESazXfffXfAWtvV23tNmoB79+7N6tWrm/KQIiIizcYYs8PXe2qCFhERaQZKwCIiIs3A7wRsjHEZY34wxrzveT3XGLPNGLPG8xjceMUUERFpXepyD/hOIAOIqrLtHmvtWw0pQElJCbt376awsLAhYVqs0NBQEhISCAkJae6iiIhIE/IrARtjEoBLgdnAb5wswO7du4mMjKR3794YY5wMHfCstRw8eJDdu3eTmJjY3MUREZEm5G8T9DPA74Dy47bPNsasM8Y8bYxp7+2DxpiZxpjVxpjVOTk5J7xfWFhIly5d2lzyBTDG0KVLlzZ79S8iEgji4uIwxpzwiIuLa9Tj1pqAjTETgP3W2u+Oe+t+oD8wHOgM3Ovt89baFGvtMGvtsK5dvQ6FapPJt0JbrruISCDIzs6u03an+HMFPAqYaIzZDiwALjDGvG6t3WfdioA5wIhGLKeIiEirUmsCttbeb61NsNb2Bq4GPrXWXmuMiQcw7ku4y4D1jVrSRnLw4EEGDx7M4MGDiYuLo3v37pWvi4uL/Yoxffp03nqrQX3RRESkjWnITFjzjTFdAQOsAW52pkg1S5ufRuqsVPJ25hHdM5pxs8eRNDWp3vG6dOnCmjVrAHjggQeIiIjg7rvvdqq4XpWVleFyuRr1GCIiEtjqNBGHtfYza+0Ez/MLrLVJ1trTrbXXWmsLGqeI/5E2P40lM5eQtyMPLOTtyGPJzCWkzU9z9DgvvfQSw4cPZ9CgQVxxxRUcPXqU/Px8EhMTKSkpAeDw4cPVXldITU1lyJAhJCUlcf3111NUVAS4p+G89957GTp0KG+++aaj5RURkfrJyMhotmMH3ExYc8+fy5q57ivSspIy5p4/l3WvrwPg3/f/m5Kj1RNeydESPr7rYwCOHjjK3PPnkrkkE4CCrPp9J5g8eTKrVq1i7dq1DBgwgFdeeYXIyEjOP/98PvjgAwAWLFjA5MmTq43fLSwsZPr06SxcuJC0tDRKS0t54YUXKt/v0qUL33//PVdffXW9yiUiIs7q0KGDz3kYYmNjG/XYAZeAa3J492Gv248ePOrocdavX895551HUlIS8+fPZ8OGDQDceOONzJkzB4A5c+YwY8aMap/LzMwkMTGRvn37AnDdddexfPnyyvevuuoqR8spIiL1c+TIEcDdOllUVIS19oRHVlZWo5ahSVdD8sf0z6ZXPneFuKq9ju4Z7W5+Pk50z2gAwmLCqu0fERdRvzJMn86iRYsYNGgQc+fO5bPPPgNg1KhRbN++nc8++4yysjJOP/30OsUNDw+vV3lERMQ5+fn5nHfeeUycOJGHHnqo2YaDtqgr4HGzxxESVr2pICQshHGzxzl6nPz8fOLj4ykpKWH+/PnV3ps2bRpTpkw54eoXoF+/fmzfvp3NmzcDMG/ePMaMGeNo2UREpGHCw8MZM2YMo0aNatZytKgEnDQ1ieSUZKJ7RYOB6F7RJKckN6gXtDcPP/wwZ511FqNGjaJ///7V3ps6dSq5ublcc801J3wuNDSUOXPm8Mtf/pKkpCSCgoK4+eYm6RwuIiK1KC4uJjc3l6CgIJ599lkuvvjiZi2PsdY22cGGDRtmV69eXW1bRkYGAwYMaLIyNNRbb73F4sWLmTdvnmMxW9rPQESkJZo2bRo//PADq1atIjQ0tEmOaYz5zlo7zNt7AXcPOJDdcccdfPTRR3z44YfNXRQREamjGTNmMGzYsCZLvrVRAq6D5557rrmLICIidZSZmUm/fv0YO3YsY8eObe7iVGpR94BFRETq4tVXX+X0009n1apVzV2UE+gKWEREWq0rrriC7OxszjzzzOYuygl0BSwiIq3OihUrKC0tJTo6mvvvv5+goMBLd4FXIhERkQbYvHkz5513HrNnz27uotRITdCAy+UiKek/Y4mvvvpq7rvvPkdir1mzhr1793LJJZc4Ek9ERGp2yimnMHfuXJKTk5u7KDVqUQk4Li6O7OzsE7bHxsY2aM7ODh06VC5J6LQ1a9awevVqJWARkUa2Zs0awsPDOfXUU5k6dWpzF6dWLaoJ2lvyrWl7Q+Tl5dGvXz8yM90rK11zzTW89NJLANxyyy0MGzaM0047jT/96U+Vn1m1ahXnnHMOgwYNYsSIEeTl5fHHP/6RhQsXMnjwYBYuXOh4OUVEBMrLy7n22muZOnUqTTnBVIN4WwGisR5nnnmmPV56enq112PGjLFz5syx1lpbXFxsx4wZY+fNm2et+yfq82GttTk5OXbMmDH2vffes9Zau2/fvhOO501QUJAdNGhQ5WPBggXWWmv/9a9/2bPPPtv+4x//sBdffHHl/gcPHrTWWltaWmrHjBlj165da4uKimxiYqL99ttvrbXW5uXl2ZKSEjtnzhx722231Xj8438GIiJSdxkZGTYzM7O5i1ENsNr6yIl+N0EbY1zAamCPtXaCMSYRWAB0Ab4DfmWtLXb260HT8NUEPX78eN58801uu+021q5dW7n9n//8JykpKZSWlrJv3z7S09MxxhAfH8/w4cMBiIqKarLyi4i0JY1xOzJtfhqps1LJ25lHdM9oxs0e5/g6A8eryz3gO4EMoCKzPA48ba1dYIz5G3AD8IKvD/urYuk/gJCQkGqvaxMTE1Nt/7i4uAaVpby8nIyMDMLCwsjNzSUhIYFt27bxl7/8hVWrVtGpUyemT59OYWFhg44jIiL+c/p2ZNr8NJbMXELJ0RIA8nbksWTmEoBGTcJ+3QM2xiQAlwIve14b4ALgLc8urwGXNUYBm9PTTz/NgAEDeOONN5gxYwYlJSUcPnyY8PBwoqOjyc7O5qOPPgLcSxHu27evcraV/Px8SktLiYyMJD8/vzmrISIiNUidlVqZfCuUHC0hdVZqox7X305YzwC/A8o9r7sAh6y1pZ7Xu4HuDpftBLGxsXXa7q9jx44xePDgysd9991HZmYmL7/8Mk8++STnnXceo0eP5pFHHmHQoEEMGTKE/v37M2XKlMr1JNu1a8fChQu54447GDRoEOPHj6ewsJCxY8eSnp6uTlgiIgEqb2denbY7pdYmaGPMBGC/tfY7Y8z5dT2AMWYmMBOgZ8+edS5gVQ0ZalSTsrIyr9szMjIqnz/11FOVz+fOnet1/+HDh7NixYoTtgfiHKQiIuIWGR9J/t4TWyqje0Y36nH9uQIeBUw0xmzH3enqAuBZoKMxpiKBJwB7vH3YWptirR1mrR3WtWtXB4osIiJtkW2k4UXjnxiPq52r2raQsBDGzR7XKMerUGsCttbeb61NsNb2Bq4GPrXWTgWWAb/w7HYdsLjRSikiIm1aVlYWZ599Nl988YVjtyM3vLmBg5sOkjQ1iUmvTiK6VzQYiO4VTXJKckD1gj7evcACY8wjwA/AK/UNZK3F3a+r7Wmsb3QiIq2JMYaysjLKysocuR1ZlF/ER3d8RJ/xfbh83uUkTU1q9IR7PNOUCWDYsGF29erV1bZt27aNyMhIunTp0uaSsLWWgwcPkp+fT2JiYnMXR0Qk4OTk5BATE4MxhvLyckdXNTrw4wGie0UT0iHEsZjHM8Z8Z60d5u29Zp8LOiEhgd27d5OTk9PcRWkWoaGhJCQkNHcxREQCzoEDBxg2bBhTp07l0UcfdST5fvPUN7jauxhx2whi+sc4UMr6a/YEHBISoqs/ERE5QZcuXZg+fTqXXebMNBO23LLj8x0Edwhm+K3Dm73VtdmboEVERKrasGED0dHRjrUOWmspKyojODSY0qJSglxBBAU3zVpENTVBt6jVkEREpHUrKSlh4sSJTJs2zbGY/77338z/+XxKC0sJbh/cZMm3Ns3eBC0iIlIhJCSE119/HSfnjYgbHEdZSRmu9q7ad25CaoIWEZFm9/XXX7Nz506uvvpqR+KVl5Xz06afmr2jlZqgRUQkoP3P//wPDz30ECUlJbXv7IdPZ33KSyNe4vCew47EawxqghYRkWb3xhtvUFBQQEiIM2NyR9wxguhe0UR1D9y12XUFLCIizeKDDz7gqquuori4mIiIiAav4V5aWMr3r3yPtZao7lEMv2W4QyVtHErAIiLSLHbt2sXWrVs5cuSII/HWzF3DkhuXsOdbr2sDBRx1whIRkSZVUFBAREQE4B525FSzs7WW3d/spsc5PRyJ54SAnopSRERar7i4OLKzs0/YHhMTQ05OTr2Tb9r8NFJnpZK3M4+QsBDGPzGe4bcOD6jkWxslYBERaTTeki+453mur7T5aSyZuYSSo+4e0yVHSvjkN58QGh3a5CsaNYTuAYuISIuSOiu1MvlWKCsqI3VWajOVqH6UgEVEpEXJ25lXp+2BSglYRERalKgE72N7o3tGN3FJGkYJWEREHLd//36mT5/uaMzMJZmUFZdx4WMXEhJWvfNWSFgI42aPc/R4ja3WBGyMCTXGfGuMWWuM2WCMedCzfa4xZpsxZo3nMbjxiysiIi1BRkYGixYtolOnTl7fj42NrVO8rLVZLJi4gJXPrSRpahLJKclE94oGA9G9oklOSW5RHbDAj3HAxr1icbi1tsAYEwJ8CdwJ3Ay8b619y9+DaRywiEjrlpubW5l0Dx06RMeOHR2LvfGDjfS5qA+ukMBa1agmDVqMwboVeF6GeB5NN3uHiIi0CJ9++im9e/dm+fLlAA1OvoWHCll4+UL2b9gPQN9L+7ao5Fsbv+4BG2Ncxpg1wH5gqbV2peet2caYdcaYp40x7X18dqYxZrUxZnVOTo5DxRYRkUAzdOhQJk+ezMCBAx2JV3S4iH0/7ONARv3HDAeyOk1FaYzpCLwL3AEcBLKAdkAKsMVa+1BNn1cTtIhI61JWVsZrr73Gddddh8vlzNVp3q48ohKiMMZQWlhKcGjLnTPKsfWArbWHgGXAz6y1+zzN00XAHGBEw4sqIiItyccff8wNN9zA4sWLHYmXk57D/w34P7578TuAFp18a+NPL+iunitfjDEdgPHAj8aYeM82A1wGrG/MgoqISOCoaD299NJL+fzzz5k8ebIjcWP6x3D2XWfTb2I/R+IFMn+ugOOBZcaYdcAq3PeA3wfmG2PSgDQgBnik8YopIiKBYsuWLYwePZqtW7cCMHr06AbFs+WWr5/8mqMHj2KCDBc8cgGRJ0U6UdSAVuu1vbV2HTDEy/YLGqVEIiIS0EpKSsjOziYnJ4eTTz65wfEOZB7g01mfYoxh5G9GOlDClkHrAYuIiF82btxI3759ASgtLSU4uGH3Z8vLyglyuRticzJyiOkfg/uuZuvhWCcsERFpmz799FMGDBjAO++8A9Dg5Ju3K48Xh7zI5k82A9B1QNdWl3xr03q7l4mISJ3FxcV5XcM3NjaWBx98kIsvvtiR44RGhxLWJaxV93KujZqgRUSkUk1XofXNF2nz00idlUrezjzCu4Uz/onxDJo2CGttq7/qrakJuu1+9RARkUaXNj+NJTOXUHK0BIAj2Ud478b3CHIFtbjFE5yme8AiItJoUmelVibfCuUl5aTOSm2mEgUOJWAREWk0eTvz6rS9LVECFhGRRuNrQo3ontFNXJLAowQsItLGHThwgIcffpjy8nJiY2O97uNruzfWWrZ/th2A8Y+PJzisenejkLAQxs0eV+/ythZKwCIibdzixYuZPXs2a9euJSsrC2vtCY+srCy/42W8ncFrY19j04ebSJqaxMSUiUT3igYD0b2iSU5JbvMdsEDDkERE2qwDBw4QExODtZatW7fSp0+fBsWrGFZUXlbOunnrGDRtECaodQ8zqo1mwhIRkWoef/xxkpKSyMrKwhjT4OS7bdk25pw3h8K8QoJcQQyePrjNJ9/aaBywiEgblJyczE8//URMTIwj8YJcQRQXFFOYW0hodKgjMVs7NUGLiLQRK1eu5Msvv+S3v/2tI/GKC4rZ+dVOTrn4FMC9rKCueqtTE7SIiPDaa6/x/PPPU1BQ4Ei81FmpLLxsIQVZ7nhKvnWjK2ARkVasqKiIQ4cOERsbS1FREQUFBXTp0qVBMSuWESw8VEjW2ix6j+ntTGFbIV0Bi4i0QdZaJk2axIQJEygrK6N9+/YNTr7LH1nOgkkLKC8rJ7RjqJJvA9TaCcsYEwosB9p79n/LWvsnY0wisADoAnwH/MpaW9yYhRUREf8ZY7jtttsoKirC5XI5EjOsaxjhXcMpL3VfBUv91doEbdxrRYVbawuMMSHAl8CdwG+Ad6y1C4wxfwPWWmtfqCmWmqBFRBqXtZbnn3+erl27cuWVVzoSMzstm6K8Inqe27NyScLWvoygUxq0HKF1/7Qr7tiHeB4WuACY4tn+GvAAUGMCFhERZ8XFxZGdnX3C9tDQ0Hon4Krr90b3iMa4DO2j2vNf3/+XOlo5yK/2A2OMyxizBtgPLAW2AIestaWeXXYD3X18dqYxZrUxZnVOTo4TZRYREQ9vyRegsLCwXvEq1u/N25EH1r1qUcG+AobeNFTJ12F+JWBrbZm1djCQAIwA+vt7AGttirV2mLV2WNeuXetZTBERaQre1u8tLSzl6z9/3Uwlar3qdAfdWnsIWAaMBDoaYyqasBOAPQ6XTUREatAYw0i1fm/TqTUBG2O6GmM6ep53AMYDGbgT8S88u10HLG6sQoqISHUbN25k6NChjsUrOlzE+ze/T2S81u9tKv5cAccDy4wx64BVwFJr7fvAvcBvjDGbcQ9FeqXxiikiIlV17tyZ4GDnpvMvPFRI+lvpnHrpqYSEhVR7T+v3Ng7NhCUi0kJkZ2eTkpLC//t//w9jDNZa4uPjvXbEio2NrXUN3/KycjZ9sIl+E/sB7iQc2jG0ei/ontGMmz1O6/fWU03DkJSARURaiL/97W/8+te/5ttvvyUpqeEJ8buXvuP9me8z48sZ9BzV04ESyvGUgEVEWqjCwkK2bdvGgAEDsNaybds2Tj755IbFzHMvGVhWUsbmjzbTN7mvJtZoJJoLWkSkhZo2bRrjx4/n2LFjGGManHz/fd+/efmslyk+UowrxEW/if2UfJuJc3fwRUTEEdZaysvLcblc3H///ezdu5cOHTo4ErvPRX0ICg7C1c6ZuaGl/tQELSISQIqLi7nqqqtISkrioYceanA8W2758vEvCe0YyvBbhjtQQqkLNUGLiLQQ7dq1IzY2tsHLBlYysOvLXexZobmSAo2ugEVEmllxcTGPPPIIN9xwA7169XIk5o+LfqTneT0J6xJGybESgkODda+3GTRoNSQREXGOr9WLjDF07tyZu+66q84xjx+3O/I3I/nX3f9i5G9HcuFjFxLSIaT2INLkdAUsItKEaroKrc//xxWrF1VdQCEkLISzf3s2Y/4wBleIOls1J90DFhFppbytXlRytIR1f1+n5BvglIBFRJqIVi+SqpSARUSayO233+5ovK+e+MrneF6tXhT4lIBFRBpRxaQaAJdeeqmjsSPiIjhpxEkEd6jen1arF7UMSsAiIo3k8OHDXHTRRbz44osAXHLJJcTGxnrd19f2qkqOlfDRf39E+lvpAAyaNojrl1/PxJcmEt0rGgxE94omOSVZqxe1ABqGJCLSSCIjI4mKiiI0NLRyW21LBNbEFeJi94rddOjcgYG/GFi5PWlqkhJuC6QrYBERB23cuJHLL7+c3NxcjDG8/fbbzJgxo97xjv10jKW/W0rJ0RKCgoOY8cUMzn/gfOcKLM1GCVhExEFHjhzhm2++IT093ZF42euyWfHMCnZ8sQOA4PZquGwtak3Axpgexphlxph0Y8wGY8ydnu0PGGP2GGPWeB6XNH5xRUQCz/fff88LL7wAwJAhQ9i+fTujRo2qd7z8vflkvpcJQO/ze3Pntjs55eJTHCmrBA5/roBLgd9aawcCZwO3GWMqbj48ba0d7Hl82GilFBEJYM8//zyPPfYYR48eBah2z7c+lv5uKYtnLKa4oBiAqO5RDS6jBJ46T0VpjFkM/C8wCiiw1v7F389qKkoRaS2++OILEhISSExMJC8vD2stHTt2rHe83K25tItoR3i3cPL35lN8pJgupzq0IpI0G8cWYzDG9AaGACtxJ+DbjTHTgNW4r5JzvXxmJjAToGfPnnUquIhIc/K1cEK3bt0oKirisssuY+7cuURH123Si+MXTxj9h9EsvWcp/Sf1Z9KcSUSeFOlUFSSA+X0FbIyJAD4HZltr3zHGxAIHAAs8DMRba6+vKYaugEWkJalp4YQvvviCwYMHExERUaeYvhZPGHLDEEb9bhRRCWpubk0avBiDMSYEeBuYb619B8Bam22tLbPWlgMvASOcKrCISKA799xz65x8wffiCZnvZSr5tjH+9II2wCtAhrX2qSrb46vsdjmw3vniiYi0HtZaLZ4glfy5BzwK+BWQZoxZ49n2e+AaY8xg3E3Q24H/apQSioi0EounL8YV4qKsuOyE97R4QttTawK21n4JeLsRomFHItLq7N27l9jYWFwuZ9bSLcovol1EO4wxnPLzUygvL+fHd3484R6wFk9oezQTloiIR2ZmJqeeeiopKSmA7wUS/Fk4AeDgxoM8d+pzbPjnBgBOv/p0Js+bTHJKshZPEC3GICJtm7WWffv2cdJJJ9G3b1/uuecefvaznwH1XzihuKCYdhHt6NSnE32T+54wnleLJwjoClhE2rhZs2YxePBgDh06hDGGBx54gMTExHrHW/7Icl444wVKjpUQ5Api4ksTiR8aX/sHpc3RFbCItDn5+fkYY4iIiOCXv/wlMTExhIeH1zteWUkZWHC1c9Hz3J4UHS7CltdtlkFpe+o8FWVDaCIOEWluhw8fZuDAgUyZMoUnnniiwfEKDxXyyjmvcMavzuC8+89zoITSmjR4Ig4RkZYuJycHgKioKG6//XYmT57coHgVvZhDO4bS56I+xA2Oa3AZpW3RFbCItBq+5m6OioqipKSEtLQ0+vTpU6eYx8/bPG72OAiCT+76hJvX3kxEXN1nw5K2Q1fAItImeEu+4G52vuWWW+jSpW6rC1XM25y3Iw8s5O3IY8nMJeTvzafPRX28z5Ag4iddAYtIq1HT4gn1+b/umd7PuJPvcaJ7RXPX9rvqHE/aHl0Bi4jUg+ZtlsakBCwircKiRYscjZf+drp7pnsvNG+zOEEJWERarKysLNLT0wEYP358g+Md++kYh7YfAiBxbCJ9ft6H4A7Vp0vQvM3iFCVgEWmRrLWMHTuWmTNnAhAeHt6guZttueWlES/x/s3vA9Chcweu/fBaJr40UfM2S6NQJywRaTGOHTvGG2+8wYwZMwgKCmLZsmX06NGDU045pV7xykrKyFycyYArBmCM4cfFP9Lp5E7EJvm32IJIbdQJS0RahSVLlnDjjTeybNkyAMaOHVvv5Auwbt463vzlm+z6ahcA/Sf1V/KVJqO5oEUkYFlreeeddwgODmbSpEn84he/4KuvvuKcc86pd7zNH28muH0wiRckcsa1ZxDZPZIeo3o4XHKR2ikBi0hAe/zxx+ncuTOTJk0iKCio3skX3Pd5//Wbf9GpTycSL0jE1c7FKRfX/wpapCFqTcDGmB7A34FY3J3yU6y1zxpjOgMLgd7AduBKa21u4xVVRFobX1NHhoaGcuDAAcLDw1m0aJFfnagqHD915PDbhpO7JZefP/dzXCEurnn/GqJ7aBiRND9/7gGXAr+11g4EzgZuM8YMBO4DUq21pwKpntciIn7zNXVkYWEhGzZsAOCkk07C5XL5Fc/b1JHL/rCMNX9fw/71+wHo3Kczrnb+xRNpTLUmYGvtPmvt957n+UAG0B2YBLzm2e014LLGKqSItD0jRoyo82dSZ6VWrlJUoayojPAu4cQPiXeqaCKOqFMvaGNMb2AIsBKItdbu87yVhbuJ2ttnZhpjVhtjVlcsByYicuDAAUfjFeYV+pwi8vCew44eS8QJfidgY0wE8DZwl7W22m+zdQ8m9jqg2FqbYq0dZq0d1rVr1wYVVkRah5SUFHr16uVYvJV/Xcmzic8S1T3K6/uaOlICkV8J2BgTgjv5zrfWvuPZnG2Mife8Hw/sb5wiikhrsGLFCrZv3w7AuHHjuPPOOxsUb/fK3RRkFQDQa0wvzvyvMxn9x9GEhIVU209TR0qgqjUBG/f6Xq8AGdbap6q89R5wnef5dcBi54snIq1Bbm4uY8eO5c9//jMAffr04dFHH6331JH5e/N5ddSrrHh2BQBxg+K48LELOfOmM0lOSdbUkdIi1DoVpTHmXOALIA0o92z+Pe77wP8EegI7cA9D+qmmWJqKUqTtWLp0KcuXL+fhhx8G4NNPP2XEiBFERETUK97mTzaTtSaLc+89F4CNH2yk95jetIto51iZRZzWoKkorbVfWmuNtfYMa+1gz+NDa+1Ba+04a+2p1toLa0u+ItL6lZeXVz7/8ssvef3118nPzwfgggsuqHPyrXqBsPmjzfzw8g+UFpUC0PfSvkq+0qJpLmgRcURGRgannXYaX375JQD33nsvGzduJDIysl7xstZk8fzA58lOc48VHvvQWG5Nv5Xg9prAT1oH/SaLiF98zVoVExNDTk4OvXr1onv37pSVlQEQFhbmV9yqM1dFJURxzu/O4azbzyK6ZzRhMWEU5xcD0D6qvXOVEQkAWo5QRPzi7o/pXXl5eY3v+1Ixc1XVyTNMkOHyv1+ujlPSKmg5QhEJSEvvXXrCzFW23JI6K7WZSiTSdJSARaRWfoyWqFMsW+6Ol78n3+s+vma0EmlNlIBFxKfDhw8zevRoXnjhBUfiFWQV8PKIl0n7RxoAUTsxGdMAABiLSURBVD00c5W0XUrAIlLNoUOH+OqrrwCIjIwkLi6u3j2ZAY7sP8Kub3YBEN4tnIj4iMrZqi587ELNXCVtljphiUg1V155JcuWLWP37t20b/+fnse+ekHHxsaSlZXlM968i+ZxcONB7tx6JyboxKbq49fvHTd7nDpgSatRUycsJWCRNm79+vXMmjWLl156iW7durF+/XpKSkoYMmRIveLtWbWHzx/8nCv+cQXtI9uz74d9hHQIIaZ/jMMlFwl8NSVgjQMWaYOOHj3K0aNHiYmJITg4mG+//ZYff/yRbt26cfrpp9c5XmFeIbbc0qFTB7CwP20/P236ifih8VqHV8QHXQGLtFI1NRm7XC4mTZrE888/D0BpaSnBwbV/H/fWXNw3uS9P93yaYbcM48LHLgSgvKycIJe6mIjoClikDfKWfCu2P//885xxxhmV2/xNvlUnzcjbkceSmUtITklm7ENj6Xlez8p9lXxFaqcrYJFWqqaxufX5u3+6x9Mc3n34hO3RvaK5a/tddY4n0hZoJiyRNmb58uWOxvtx0Y9eky9o0gyR+lICFmkFiouLefbZZ/nkk08AGDp0aMPiFRTz+UOfs/3z7QD0GNWDdpHel/7TpBki9aMELNJCWWvZs2cP4L6H+8wzz7BkyRKAei16X1ZSxqEdhwBwtXOx6vlV7Ph8BwDhXcOZ8MIETZoh4qBae14YY14FJgD7rbWne7Y9ANwE5Hh2+7219sPGKqSInGjatGl88803bNy4kaCgIL777js6d+5c+X5sbKzPXtDevHHpGxz76RgzV8/E1c7FHRvvqLYEYMXkGJo0Q8QZtXbCMsaMBgqAvx+XgAustX+py8HUCUuk/tauXcuf//xnXnzxRcLDw1m6dCn79u1jypQpfvViPt6OL3bw7V+/ZfIbk3GFuNj04SbKy8rpO6FvvZYWFJETNWgYkrV2uTGmt9OFEpHqfI3b7dq1K/v37+fw4cN8/PHHpKenM3z4cMaPH19rzOPH7Q6/dThnzjyT0I6hFB4qZO/qvRzadogufbtw6iWnNka1RMQHv4YheRLw+8ddAU8HDgOrgd9aa3Nri6MrYBHfahs2ZK2luLi42vzMNfG22D24m5Invz65cklAb/Mzi4gzGmMY0gtAH2AwsA94soaDzzTGrDbGrM7JyfG1m4jUwhjjd/IFSP196gnJF6jsWGWCjJKvSDOqVwK21mZba8usteXAS8CIGvZNsdYOs9YO69q1a33LKdIqbd++nVdeecWxePt+2Mfqv7lbmfJ2eR+fe3iP9/G8ItK06pWAjTFVZ1e/HFjvTHFEWr8DBw5QVlYGwIIFC7jpppvYt29fvePlbs2tnNlq/YL1LL1nKcVHin2Oz9W4XZHAUGsCNsb8A/gG6GeM2W2MuQF4whiTZoxZB4wFft3I5RRpFZYvX058fDyfffYZADfddBPbt28nPr5+Kwalv5XOX/v8lawf3OvxjrpnFL/e/Wvahbdj3OxxGrcrEsBqTcDW2mustfHW2hBrbYK19hVr7a+stUnW2jOstROttfX/+i7Sih05coRp06bx+uuvAzB8+HDuueceTj75ZAC6dOlCz57uRQx8jc+tuv3I/iO8/rPXSX87HYDEcYmM//N4onpEARAWE0ZodCjg7myVnJJMdK9oMO45m5NTkjVuVyRAaDUkkXryNWyoc+fOzJ07l+TkZMLCwti0aRODBw8GoEOHDjz66KNe42VlZXld7i+0UyiZSzLpl9yPDl06UJxfTGlhqTtepw6cc/c5PsuYNDVJCVckQGk1JJF6qmnY0CmnnMLGjRsxxmCt9WtiC2/DhkLCQojsHklEXAQzls9wpNwi0nS0GpJIE/v2228rk64/yddaS+qsE4cNlRwtoeRoCdP+Pa1RyikizUcJWKSOvv/++1pXG+rUqZPf8TZ/spmnE572uaxf/t58XO1cdSqjiAQ+JWCRWhQUFPCHP/yB1NRUAE466SRCQ0PrHy+rgMU3LGbX17sA6JTYiYSzE4iI976CkYYNibROSsAiXqxcuZJly5YBEBoaSkpKCitWrADcna++/vprv2PZckvGOxls+3QbAO0i27H5w80c3HQQgC59u3Dl21dy0RMXadiQSBuiXtAiQFlZGVu2bKFv374A3HHHHQQHB/P1118THBzMtm3bCAsLq/aZmpb7K8wrJHdrLvFD4sHA0t8tJW5wHIkXJNIuvB2/3v1rglzVv/9quT+RtkW9oKVN8DVkKDY2lqysLGbOnMk777xDVlYWwcHBpKenc9JJJ9GxY0efMX31Wk5OSWb9gvXkpOdwx+Y7MMaQuy2X6B7RBAWr0UmkLWnQcoQirYG35Ft1+w033MBFF11UOaXjwIEDa43pq9dy6qxUfrHgF5WrDYH7Pq+ISFVKwCLAWWedxVlnneX3/rtX7iZvh/dey3k780g4O8GpoolIK6X2MGmVjh07xj333MOiRYsciXck5whvXvkmmz7aBEBU9yhc7b0PDVKvZRHxhxKwtBrz58/njTfeANw9l999913Wrl1br1jWWj5/+HPWvb7OHa9jKPvT9nMk+wgAUQlRTHplknoti0i9qQlaAlJtnaYA1q5dS3p6Otdccw0Ac+bMwVrLlClTMMbw448/Ehzs/6/4pg83UZBdwJAZQzDGkLk4k5OGn8QZ156BK8TFbRm3VdtfvZZFpCHUC1oCUk3TN1b8zt56663Mnz+fAwcOEBISwk8//USnTp28fjamYwwH8w6esL1jh47kHs0F4K2r3yJ7bXZloi0rLtMMVCLSIDX1glYCloBUUwLesmULJ598Mnv27KF9+/bExMTUGu+Z3s/47DR1b+69hHYM5eiBo4R2DNVQIRFxjIYhSYuxc+dOwsPDa9ynYv3c7t2717hfWXEZu1fupttp3XzOs4xx398F91q6IiJNRV/1pVnl5uayd+9eAHbv3k2vXr0qF6/3xdd9XVtuyV6XXZls96/fz9zRc9n04SafPZPVY1lEmkutCdgY86oxZr8xZn2VbZ2NMUuNMZs8/2qWAfFLYWEhO3fuBNzTPyYmJvLII48AkJCQwIsvvsiECRP8jndoxyEO/HgAgOKCYl4c+iLfvfQdAHGD47jq3avoO6Ev42aPU49lEQko/lwBzwV+dty2+4BUa+2pQKrntbRhcXFxGGNOeMTFxbFjx47K/UaPHs31118PgMvl4n//938rXwPMnDmTPn36EBsb6/U43bp2I3udu3e0tZY5585h2R/diya0j2rPVe9exfBbhgNgggz9L+tPaMdQkqYmkZySTHSvaDAQ3Sua5JRk9VgWkWbjVycsY0xv4H1r7eme15nA+dbafcaYeOAza22/2uKoE1brVVOnqYSEBHbu3IkxhsWLFxMaGsrFF19cYzxf8yxH947GFezi5rU3A7Bl6Raie0QT07/2jlgiIk2tMTphxVpr93meZwHeL1fcB58JzIT/dJ6R1sNaS0ZGRo37PP7445SXl+NyuZg0aVKN+xYeKmT3it2k/t77PMvHDh7jqnevqtzWZ3yf+hdeRKQZNbgTlnVfQvu8jLbWplhrh1lrh3Xt2rWhh5NmVlpayqpVq8jNdY+d/ec//8lpp51W42emTJmCy+V9PO2xn46xfuF6iguKAVg7by3zfz6fvF3eey0f2X+EHiN7NKAGIiKBob4JONvT9Izn3/3OFUkCSUlJCcuWLWPLli2Ae/apESNG8PHHHwMwduxYXn31Vb/jFR4q5PtXvid3mzuB7/1uL29f/Ta7vt4FwMArBnLdZ9cR1SPK6+fVa1lEWov6JuD3gOs8z68DFjtTHGkKNXWYKi8v54MPPuCbb74B4OjRo4wbN4558+YBMGjQIBYuXMj48eMB6NatGzNmzKjxeCufW1mZYAvzClly4xI2f7wZgJ6jenLTqptIvCARgMiTIuk9pjcXPnqhei2LSKtWaycsY8w/gPOBGCAb+BOwCPgn0BPYAVxprf2ptoOpE1ZgqKnDVHl5OQkJCYwZM6ZyYYPPPvuMQYMG0amT79FmUa4o8svzT9geGRTJfR3u46z/Potxj7qT54HMA3Tp26XGcoC7I5bmWRaRlkxTUUo1tc2znJ6eTu/evQkL829mqM8e/IzPH/jcx8Hg7uy7Ce9a8+xWIiKtUU0JWDNhtVJVv1jNmTOHiRMn+v3ZgQMHnpB8i48UVz5PnZXKvPHzKl/n7cwjJLx6c3GF6J7RSr4iIl4oAbcANd2zrZCbm0tZWRkA8+bNIy4ujiNH3GvXFhUVUVBQwLFjx/w6XllJGVlrsipfL/vjMp466SnKy8oB933aTqd0qkzyk16ZRPKLybpnKyJSB0rALYC3dXGrbn/33Xfp3Lkz6enpAPTu3Zvk5GQKCgoAuPnmm/n000/p0KGDX8db9fwqXhzyIof3HHbHG9ubUfeNoqzYneBH3DaCCS9MqNaUrZmmRETqRveAW4Da7tnu2rWLefPmMW3aNBISEmqN56vDVERQBPll+eRuzWXv6r2cesmptIto16Cyi4i0ZVqOsIUpKiri7rvvZuzYsUyePLnW/Xv06MHvf//7atustZQVlxHcPpiC7AIWXbeI4bcNp19yP35rf+s9kOe7WKeTO9HpZK2vISLSmNQE3Qj8uWdbWFjIvn37Kl+PHTuW3/3udwC0a9eODz/8kMzMTL+OZ63l4MaDHNx4EIDSwlL+HPNnvnriKwA6dOrA0QNHKT1WCkB0Dy3NJyLS3HQF3Ahqu2cLMHLkSOLj4/nwww8BSEpKIjHRPRmFMYbNmzfXOk62gjGGuefP5eQLT+byv19OcGgww24ZRsJZ7uZoVzsXM1fPrNx/3KPjvC50oA5TIiJNR/eAG0Ft92wB3nrrLcLCwrjkkkt87lt8pJh24e1837Mlgnzr3r7lX1uI6hFF1wH+zbetSS5ERBqfJuJoBAUFBaxfv56zzjoLYwzPPvssTz75JNu3b/e58ABUH59b1U9bfiJ7XTYDLh8AwPs3v8/mjzdz1/a7eDDoQe/LXRj4U/mfnKiOiIg0Ak3EUQt/7tlmZGTwhz/8gbw89yo9r776KiNHjiQryz1etk+fPlxyySUcPXq0xmNVJOCt/97K21PeprzUPbZ2zdw1vPnLNyktdN+n7TuhL2f991lYa33em9U9WxGRlksJmJrv2f7www8AbN26lccee4yNGzcCkJyczHvvvUdUlHvVngkTJvC3v/2NiIiIGo91eJd7bG3+vnx2r9hNQbZ7rO6ZM8/klrRbcLVzXz33ndCXkb8ZiTGGcbPHaZILEZFWps12wjp8+DApKSlccMEFNe5XWuq+Ih0/fjwFBQWEhoYCkJiYSM/uPbHl7iva3K25fP7g55z967MJJ5wjHDkhVjj/mZLxjGvPYNCvBlW+9tUzGai8N6t7tiIirUeLS8BxcXFer1hjY2Mrm4PB3dS7YcMGQkNDOeWUUzh27BgjR47kxhtv5Pbbb8cYwz333MOTTz5Z4/GGDx8OgCk1pP8jnfih8cQPiSd3ay7P9X2OSXMmuROpcTcrD7xyII/0eoS8HScuKB/dK7qy2djfHs4VkqYmKeGKiLQiLa4Juqbm4oceeqja4vDnnnsuTz31FAAdOnSgX99+xHSOASAyMpLdm3dz6/W31ni8lX9dCYAJMiy5aQmZ77nH5kb1iOLc+86l2+ndAOiU2Inf7PkNfS/tqyZjERGpVYu7Aq7J35/6O4kxiVx//fUYY5geO53ErYmV749OG02wCYYp7tdvX/Q2Pc7pUWPMdpHuqRiDQ4O5a/tdRHaPBMAV4uKCR7w3X6vJWEREatOqEvCce+dUm7v46luvJrRjaOXrs399NuHd/nMf9vwHzyciLoLw133fsx0yY0jl67r0OlaTsYiI1KTFjQP2Z5KLunqm9zM+79netf2uesUUERFptHHAxpjtxpg0Y8waY0yLnWFD92xFRKSpOdEEPdZae8CBOH6JjY312Qu6vnTPVkREmlqLuwdcdaiRk3TPVkREmlJDhyFZ4F/GmO+MMTO97WCMmWmMWW2MWZ2Tk9PAw4mIiLQODU3A51prhwI/B24zxow+fgdrbYq1dpi1dljXrv6t1CMiItLaNSgBW2v3eP7dD7wLjHCiUCIiIq1dvROwMSbcGBNZ8Ry4CFjvVMFERERas4Z0wooF3vWMyw0G3rDWfuxIqURERFq5Jp2IwxiTA+xwMGQM0GRDoBqZ6hJ4Wks9QHUJVK2lLq2lHuB8XXpZa712gGrSBOw0Y8xqXzOMtDSqS+BpLfUA1SVQtZa6tJZ6QNPWpcWthiQiItIaKAGLiIg0g5aegFOauwAOUl0CT2upB6gugaq11KW11AOasC4t+h6wiIhIS9XSr4BFRERaJCVgERGRZtCsCdgY08MYs8wYk26M2WCMudOzvbMxZqkxZpPn306e7cYY81djzGZjzDpjzNAqsa7z7L/JGHOdj+N5jRso9TDGDDbGfOOJsc4Yc5WP4003xuR41mFeY4y50Yl6OFkXz3tlVcr4no/jtTfGLPR8fqUxpneg1cUYM7ZKPdYYYwqNMZd5OV4gnZf+nt+lImPM3cfF+pkxJtNTz/t8HK9RzotT9fAVx8vxzjfG5FU5J390oh5O1sXzXq1rq9f0txYodTHG9Dvub+WwMeYuL8drlPNSj3pM9fws04wxXxtjBlWJ1fh/J9baZnsA8cBQz/NIYCMwEHgCuM+z/T7gcc/zS4CPAAOcDaz0bO8MbPX828nzvJOX43mNG0D16Auc6nl+ErAP6OjleNOB/w3kc+J5r8CP490K/M3z/GpgYSDWpUrMzsBPQFiAn5duwHBgNnB3lTguYAtwMtAOWAsMbKrz4mA9vMbxcrzzgfcD+Zx43tsOxNRyvFp/PwOhLsf9rmXhnoiiSc5LPepxDp5cgXtRoZVVyt7ofyeO/1I28Ie3GBgPZALxVX6gmZ7nLwLXVNk/0/P+NcCLVbZX2+/4/Y+PGyj18BJnLZ6EfNz26TTSf/RO1gX/EvAnwEjP82DcM9CYQKtLlW0zgfk+4gfMeamy3wNUT1wjgU+qvL4fuL+5zkt96+Erjpft59NICdjJuuBfAvbr/43mrkuV9y4CvvLxXpOcF3/r4dneCdjjed4kfycBcw/Yc+k+BFgJxFpr93neysI97zRAd2BXlY/t9mzztf14vuI6poH1qBpnBO5vXlt8HOoKT9PJW8aYHs6UvjoH6hJq3GtBrzBemmyP/7y1thTIA7o4VYcKTp0X3N9y/1HDoQLlvPji799Ko5+XBtbDVxxvRhpj1hpjPjLGnFbf8tahDPWpS61rq+P/uWsQp84Ltf+tNOp5qUc9bsDdwgBN9HcSEAnYGBMBvA3cZa09XPU96/5q4fhYqcaI61Q9jDHxwDxghrW23MsuS4De1tozgKXAaw0quPcyOFGXXtY9pdsU4BljTB+ny+kPh89LEu5vvd60lPPS7Bw8Jz7jeHyP+/dwEPAcsKhBBa9jGepQl1rXVm8KDp6XdsBE4E0fuzTqealrPYwxY3En4HudLEdtmj0BG2NCcP+g5ltr3/Fszvb8Z1fxn95+z/Y9QNWrigTPNl/bj+crbqDUA2NMFPABMMtau8Lbsay1B621RZ6XLwNnOlUPJ+ti/7Ne9FbgM9zfRo9X+XljTDAQDRwMtLp4XAm8a60t8XasADsvvvj7t9Jo58WheviKU4219rC1tsDz/EMgxBgT40A1aipDneti/Vtb3d9zVy9O1cXj58D31tpsb2825nmpaz2MMWfg/nudZK2t+B1vkr+T5u4FbYBXgAxr7VNV3noPuM7z/Drc7fgV26cZt7OBPE+zwifARcaYTp7ebRfh/SrFV9yAqIfnW+O7wN+ttW/VcLz4Ki8nAhlO1MMT26m6dDLGtPfEjAFGAeleDlk17i+ATz3fUAOmLlU+dw01NKkF2HnxZRVwqjEm0fP7drUnxvEa5bw4VY8a4hy/X5xn34rbOkE490XCqbr4u7Z6bb+f9ebg71eF2v5WGuW81LUexpiewDvAr6y1G6vs3zR/J/7eLG6MB3Au7qaAdcAaz+MS3G3oqcAm4N9AZ8/+Bvg/3PdF04BhVWJdD2z2PGZU2f5yxX6+4gZKPYBrgZIqMdYAgz3vPQRM9Dx/DNiAu5PWMqB/oJ0T3L0L0zxlTANuqHKMqnUJxd1MtRn4Fjg50Oriea837m+7QccdI1DPSxzu+1aHgUOe51Ge9y7B3Tt0C+6WliY7L07Vw1ccz2duBm72PL+9yjlZAZwTaOcEd0/btZ7HhuPOSdW6+Pz9DJS6eN4Lx51Mo487RqOfl3rU42Ugt8q+q6vEavS/E01FKSIi0gya/R6wiIhIW6QELCIi0gyUgEVERJqBErCIiEgzUAIWERFpBkrAIiIizUAJWEREpBn8f01MVrOu+zn6AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "## PLOTTING METHOD\n", "y=6*np.exp(0.1*(t-2000)) # EXACT SOLUTION\n", "fig = plt.figure(figsize=(8,4))\n", "plt.plot(t,w,'o:',color='purple',label='Taylor')\n", "plt.plot(t,y,'s:',color='black',label='Exact')\n", "plt.legend(loc='best')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": { "id": "dTcc7pw5aNbx" }, "source": [ "## Table\n", "The table below shows the time, the numerical approximation, $w$, the exact solution, $y$, and the exact error $|y(t_i)-w_i|$ for the linear population equation:" ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 696 }, "id": "R9kr9RdCaNbx", "outputId": "e76c795a-475a-411d-a7b2-9b9933cb36d1" }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
time t_iAdams wExact y(t_i)Exact Error |w_i-y_i|
02000.06.0000006.0000000.000000
12001.06.0600006.6310260.571026
22002.06.6690007.3284170.659417
32003.07.3663508.0991530.732803
42004.08.1378528.9509480.813096
52005.08.9902139.8923280.902115
62006.09.93185210.9327131.000861
72007.010.97211912.0825161.110397
82008.012.12134513.3532461.231901
92009.013.39094014.7576191.366678
102010.014.79351416.3096911.516177
112011.016.34299418.0249961.682002
122012.018.05476819.9207021.865934
132013.019.94583322.0157802.069947
142014.022.03497024.3312002.296230
152015.024.34292426.8901342.547211
162016.026.89261429.7181952.825581
172017.029.70936032.8436843.134325
182018.032.82113336.2978853.476752
192019.036.25883540.1153673.856532
202020.040.05660344.3343374.277733
\n", "
" ], "text/plain": [ " time t_i Adams w Exact y(t_i) Exact Error |w_i-y_i|\n", "0 2000.0 6.000000 6.000000 0.000000\n", "1 2001.0 6.060000 6.631026 0.571026\n", "2 2002.0 6.669000 7.328417 0.659417\n", "3 2003.0 7.366350 8.099153 0.732803\n", "4 2004.0 8.137852 8.950948 0.813096\n", "5 2005.0 8.990213 9.892328 0.902115\n", "6 2006.0 9.931852 10.932713 1.000861\n", "7 2007.0 10.972119 12.082516 1.110397\n", "8 2008.0 12.121345 13.353246 1.231901\n", "9 2009.0 13.390940 14.757619 1.366678\n", "10 2010.0 14.793514 16.309691 1.516177\n", "11 2011.0 16.342994 18.024996 1.682002\n", "12 2012.0 18.054768 19.920702 1.865934\n", "13 2013.0 19.945833 22.015780 2.069947\n", "14 2014.0 22.034970 24.331200 2.296230\n", "15 2015.0 24.342924 26.890134 2.547211\n", "16 2016.0 26.892614 29.718195 2.825581\n", "17 2017.0 29.709360 32.843684 3.134325\n", "18 2018.0 32.821133 36.297885 3.476752\n", "19 2019.0 36.258835 40.115367 3.856532\n", "20 2020.0 40.056603 44.334337 4.277733" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "\n", "\n", "d = {'time t_i': t, 'Adams w': w,'Exact y(t_i)':y,'Exact Error |w_i-y_i|':np.abs(y-w)}\n", "df = pd.DataFrame(data=d)\n", "df" ] }, { "cell_type": "markdown", "metadata": { "id": "Q2u4famhaNbz" }, "source": [ "## 2. Non-Linear Population Equation \n", "\\begin{equation} y^{'}=0.2y-0.01y^2, \\ \\ (2000 \\leq t \\leq 2020), \\end{equation}\n", "with the initial condition,\n", "\\begin{equation}y(2000)=6.\\end{equation}\n", "## Specific 2 step Adams-Bashforth method for the Non-Linear Population Equation\n", "The specific Adams-Bashforth difference equation for the non-linear population equations is:\n", "\n", "\\begin{equation}w_{i+1}=w_{i}+\\frac{0.1}{2}\\big[3(0.2w_{i}-0.01w_{i}^2)-\n", "(0.2w_{i-1}-0.01w_{i-1}^2)\\big], \\end{equation}\n", "\n", "for $i=1,...,199$, where $w_i$ is the numerical approximation of $y$ at time $t_i$, with step size $h$ and the initial condition\n", "\\begin{equation}w_0=6.\\end{equation}\n", "To solve the 2 step method we need a value for $w_1$, here, we will use the approximation from the 2nd order Taylor method (see other notebook),\n", "\\begin{equation}w_1=6.084.\\end{equation}" ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "id": "FuwlIUq1aNbz" }, "outputs": [], "source": [ "def nonlinfun(t,w):\n", " ftw=0.2*w-0.01*w*w\n", " return ftw" ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "id": "gLr1q8FkaNb0" }, "outputs": [], "source": [ "### INSERT METHOD HERE\n", "w=np.zeros(N+1)\n", "w[0]=6\n", "w[1]=6.084 # FROM THE THE TAYLOR METHOD\n", "for n in range(1,N):\n", " w[n+1]=w[n]+h/2*(3*nonlinfun(t[n],w[n])-nonlinfun(t[n-1],w[n-1]))" ] }, { "cell_type": "markdown", "metadata": { "id": "4iBCF4SEaNb2" }, "source": [ "## Results\n", "The plot below shows the numerical approximation, $w$ (circles) for the non-linear population equation:" ] }, { "cell_type": "code", "execution_count": 23, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 295 }, "id": "u5yFmrlgaNb2", "outputId": "be9a4c7a-fe10-4080-c45b-82a2b31f3511" }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAe4AAAEWCAYAAACg1nQiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXxU1f3/8dcnIWENYVXWsAmUJYIYQcEKilXUoq3YuqDiSv25Fdva2rrb4tfuWLUqbrhQShWXotaqiGJRsKAIKC5oWV3YwxYgJJ/fH3MTJ2EmmSyTmUnez8djHpm56+fMncxn7rnnnmPujoiIiKSGtEQHICIiIrFT4hYREUkhStwiIiIpRIlbREQkhShxi4iIpBAlbhERkRSixC1SRWZ2n5ndmOg4EsXMppnZb2qw/k4z61mbMSUbM/vAzEYlOg6pn5S4JWWY2Soz22BmzcOmXWJmr8dhX6PMbF2kee5+mbv/urb3WR3Be1IQJMOvg6TaItFxlTCz183skvBp7t7C3T+Pw77C34uSx921vZ8I+z3gh4y7D3D31+O9b2mYlLgl1aQDP050EHXNQqL9v4519xbAECAPuKHuIks6Y4MfBiWPKxMdkEhtU+KWVPN74Gdm1irSTDMbbmb/NbP84O/wsHmvm9mvzWy+me0ws5fNrF1VAwg/wyo5Mzeznwa1AV+a2YVhyzY2sz+Y2ZrgjPg+M2sazGttZs+b2UYz2xo871Iu3slmNh/YDVRYvezu64F/AQOD9U8Nqmy3BdvqF7btVWb2SzP7MNj3I2bWJJh3gZn9p1yZ3cwOifBeRC2DmU0Gvg3cHX72G74tM8s2s8eC9Veb2Q0lP1BK4gjev61m9j8zOynW41QuzvRgO5vM7HMzuyKIo1HY+3F82PK3mNkTYa+fNLOvgs/VPDMbEEyfCIwHfh6UcXb57QWfgSlm9kXwmGJmjYN5FX5+RCJR4pZUswh4HfhZ+Rlm1gZ4AfgL0Bb4E/CCmbUNW+wc4ELgICAz0naqoQOQDXQGLgbuMbPWwbw7gD7AYOCQYJmbgnlpwCNANyAHKADKV+2eB0wEsoDVFQVhZl2Bk4H3zKwPMAOYBLQHXgRmm1lm2CrjgROBXkGM1TlTj1oGd78eeBO4soKz37sIvXc9gZHA+YSOT4lhwMdAO+B3wENmZtWI81Lgu8BhhGolzqji+v8CehP63LwLTAdw96nB898FZRwbYd3rgSMJfQYGAUMp+15X9PkROYASt6Sim4CrzKx9uemnAJ+6++Puvt/dZwAfAeFfpo+4+yfuXgD8g9CXaU0VAre5e6G7vwjsBPoGCWYicI27b3H3HcDtwFkA7r7Z3We5++5g3mRCySvcNHf/IChPYZT9P2tm24D/AG8E+zgTeMHdXwnW+wPQFBgett7d7r7W3bcE+z67qgWPsQwRmVk6offil+6+w91XAX8k9GOlxGp3f8Ddi4BHgY7AwRVs9tmghqHkcWkw/YfAlLDy/l8Vy/lwEONe4BZgkJllx7j6eEKfjw3uvhG4lbJljPj5qUp80rA0SnQAIlXl7svN7HngOmBF2KxOHHhWuprQmUyJr8Ke7wZqoyHXZnffH2G77YFmwOKwk0QjdJ0eM2sG/BkYA5ScYWWZWXqQqADWxrD/77n7q+ETzKzMe+HuxWa2lrLvRfi2VxN6/6okxjJE0w7IoOwxi3q83H138D5WdMwOeC8CnTiwvDEJfmBMBn5A6JgWh8WfH8Mmyn8uy7/X0T4/IhHpjFtS1c2Eqj/Dv+S/IFRlGy4HWF9XQZWziVDV8QB3bxU8soOGZAA/JXRmNczdWwLHBNPDq4KrO3xfmfciOPvvStn3omvY85xgHYBdhH5wlKzboYL9VFaGiuLfROhsM/yYxet4fcmB5Q1XpsyEqq9LnAOcBhxPqEq7ezA9ljLCgZ/L8PdapMqUuCUluftKYCZwddjkF4E+ZnaOmTUyszOB/sDz1d2PmTUp94j5+qq7FwMPAH82s4OC7XU2sxODRbIIJfZtwfX5m6sbZwT/AE4xs9FmlkEowe4F3gpb5goz6xLs+3pC7yfA+8AAMxscNFi7pYL9VFaGr4nSqC44I/8HMNnMssysG/AT4IlIy9fQP4Crg/K2JlRbE24JcJaZZZhZ+WvgWYTeu82Ekvvt5daNWsbADOAGM2tvocaQNxGfMkoDocQtqew2oPSebnffTKgB0k8Jfcn+HPiuu2+q5vY7E0pK4Y9eVdzGL4CVwAIz2w68yjfXL6cQuu68CVgAvFTNOA/g7h8D5xJq/LWJ0HX+se6+L2yxvwEvA58DnwG/Cdb9hNB7+yrwKaFr59FUVoY7gTOCVuF/ibD+VYTOdj8P9vM34OGYC3qg2Vb2Pu5ngukPAP8m9KPkXeDpcuvdSOjYbiV0DfpvYfMeI1S9vR74kFA5wz0E9A+uqT8bIabfEGpUuRRYFuy/2h3YiJh7dWviRCRVmdkq4JIo14PrPTPrDvwPyCh3fVkk6emMW0REJIUocYuIiKQQVZWLiIikEJ1xi4iIpJCU6IClXbt23r1790SHISIiUicWL168yd3L9w4JpEji7t69O4sWLUp0GCIiInXCzKL27qeqchERkRSixC0iIpJClLhFRERSSEpc446ksLCQdevWsWfPnkSHIimuSZMmdOnShYyMjESHIiJSqZRN3OvWrSMrK4vu3btThXEfRMpwdzZv3sy6devo0aNHosMREalUyibuPXv2KGlLjZkZbdu2ZePGjYkORURS0LLpy5hz/Rzy1+STnZPN6MmjyR2fG9d9pmziBpS0pVbocyQi1bFs+jJmT5xN4e5CAPJX5zN74myAuCZvNU4TERGphjm/mlOatEsU7i5kzvVz4rpfJe5qWrt2Lcceeyz9+/dnwIAB3HnnnTGvu2TJEl588cVajWfJkiWYGS+9FH1I5wsuuICnnnqqVvdb0b6aNWvGjh07SqdNmjQJM2PTpoqHx7799ttLn69atYqBAwdWO46ari8i9cuy6cuY0n0Kt6bdypTuU1g2fVmFy+/etJsdX4S+x7zYmXn6TBZPXQxA/tr8iOvkr4k8vbY0mMRd1YNVmUaNGvHHP/6RDz/8kAULFnDPPffw4YcfxrRuPBL3jBkzOProo5kxY0atbrcmDjnkEJ577jkAiouLee211+jcuXOl64UnbhGR2lJStZ2/Oh/8m6rtpU8sLV1m4V8WsuTRJaWv7829l9dufA0ASzP2bNtTepadnZMdcT/RpteWBpG4ox2smiTvjh07MmTIEACysrLo168f69evP2C5J598koEDBzJo0CCOOeYY9u3bx0033cTMmTMZPHgwM2fOZNeuXVx00UUMHTqUww47rDTZTZs2jdNOO41Ro0bRu3dvbr311oixuDtPPvkk06ZN45VXXim9Rc7dufLKK+nbty/HH388GzZsKF3ntttu44gjjmDgwIFMnDiRklHiRo0axTXXXENeXh79+vXjv//9L6effjq9e/fmhhtuAGDXrl2ccsopDBo0iIEDBzJz5syIcZ111lml815//XVGjBhBo0bfNKt44oknGDp0KIMHD+ZHP/oRRUVFXHfddRQUFDB48GDGjx8PQFFREZdeeikDBgzghBNOoKCgAAj9ADryyCM59NBD+f73v8/WrVsBWLx4MYMGDWLQoEHcc889sRxOEWkA5lwfuWp79qWzS19/MPMDPn3+09LXY+4cw5BLhpS+nvDaBI6cdCQAoyePJqNZ2dtIM5plMHry6HiEX6reJO5po6axZFroV1JRYRHTRk0r/RX16i9fjXiwXpoUqlbevWk300ZN4+PZHwOw86udVdr3qlWreO+99xg2bNgB82677Tb+/e9/8/777/PPf/6TzMxMbrvtNs4880yWLFnCmWeeyeTJkznuuON45513mDt3Ltdeey27du0C4J133mHWrFksXbqUJ598MmKf7W+99RY9evSgV69ejBo1ihdeeAGAZ555ho8//pgPP/yQxx57jLfeeqt0nSuvvJL//ve/LF++nIKCAp5//vnSeZmZmSxatIjLLruM0047jXvuuYfly5czbdo0Nm/ezEsvvUSnTp14//33Wb58OWPGjIn4vvTp04eNGzeydetWZsyYwVlnnVU6b8WKFcycOZP58+ezZMkS0tPTmT59OnfccQdNmzZlyZIlTJ8+HYBPP/2UK664gg8++IBWrVoxa9YsAM4//3x++9vfsnTpUnJzc0t/2Fx44YXcddddvP/++1U6jiJSP3jxN8NVv//4+7xweeg7MVoV9v49+0ufT3h9Aj948gelrwf8cABdj+oacb3c8bmMnTqW7G7ZYJDdLZuxU8fGvVV5vUncFdm+bnvE6bs3767xtnfu3Mm4ceOYMmUKLVu2PGD+iBEjuOCCC3jggQcoKiqKuI2XX36ZO+64g8GDBzNq1Cj27NnDmjVrAPjOd75D27Ztadq0Kaeffjr/+c9/Dlg/PCmeddZZpdXl8+bN4+yzzyY9PZ1OnTpx3HHHla4zd+5chg0bRm5uLq+99hoffPBB6bxTTz0VgNzcXAYMGEDHjh1p3LgxPXv2ZO3ateTm5vLKK6/wi1/8gjfffJPs7OjVQqeffjp///vfWbhwId/+9rdLp8+ZM4fFixdzxBFHMHjwYObMmcPnn38ecRs9evRg8ODBABx++OGsWrWK/Px8tm3bxsiRIwGYMGEC8+bNY9u2bWzbto1jjjkGgPPOOy9qbCKS3GK5xFmwpYDP53xeWms4/3fz+V2731FcVAzAlpVbWLdgHV7s0au2u30zPT0jvUox5o7PZdKqSdxcfDOTVk2Ke9KGFL8dLNwFr19Q+jw9I73M6+yc7FA1eTklB7FZu2Zllm/RoUVM+ywsLGTcuHGMHz+e008/PeIy9913HwsXLuSFF17g8MMPZ/HixQcs4+7MmjWLvn37lpm+cOHCA25VKv+6qKiIWbNm8dxzzzF58uTSDkXCG4WVt2fPHi6//HIWLVpE165dueWWW8r0QNe4cWMA0tLSSp+XvN6/fz99+vTh3Xff5cUXX+SGG25g9OjR3HTTTRH3deaZZ3L44YczYcIE0tK++Z3o7kyYMIH/+7//ixpn+XgA0tPTS6vKRaT+inarVf66fHZt2MUxNxxD09ZNWfa3Zfzrqn9xzbpraNm5JR0Gd2DIJUPYv2c/mc0zOfbWYzn21mOBUNV2+Dahbqq2a1uDOOOOx3UId+fiiy+mX79+/OQnP4m63GeffcawYcO47bbbaN++PWvXriUrK6tMYj3xxBO56667Sn8xvvfee6XzXnnlFbZs2UJBQQHPPvssI0aMKLP9OXPmcOihh7J27VpWrVrF6tWrGTduHM888wzHHHMMM2fOpKioiC+//JK5c+cClCbpdu3asXPnziq3NP/iiy9o1qwZ5557Ltdeey3vvvtu1GW7devG5MmTufzyy8tMHz16NE899VTpdfctW7awenVoFLuMjAwKCwsP2Fa47OxsWrduzZtvvgnA448/zsiRI2nVqhWtWrUqrZkoqW4XkdRRtK+IV38V+RLn2396m0V/XcTWz0JtWvqe1pfz55xPs7bNAOh1Qi++87vvkNk884DtJqpqu7bVmzPuipQclNrs3Wb+/Pk8/vjj5Obmllbj3n777Zx88slllrv22mv59NNPcXdGjx7NoEGDyMnJKa0a/+Uvf8mNN97IpEmTOPTQQykuLqZHjx6l15yHDh3KuHHjWLduHeeeey55eXlltj9jxgy+//3vl5k2btw47r33Xl588UVee+01+vfvT05ODkcddRQArVq14tJLL2XgwIF06NCBI444okplX7ZsGddeey1paWlkZGRw7733Vrj8j370owOm9e/fn9/85jeccMIJFBcXk5GRwT333EO3bt2YOHEihx56KEOGDGHy5MlRt/voo49y2WWXsXv3bnr27MkjjzwCwCOPPMJFF12EmXHCCSdUqWwiUrfcna2fbSWjWQZZnbLY/Olm7s29l6K9kS8t7t64mxsLbyQtPXTemd01m+yusbfizh2fm3KJujwrOctLZnl5eV6+UdaKFSvo169fgiKqG9OmTWPRokXcfffdiQ6l3msInyeReIul+093Z+W/VtK0TVO6HNmFvTv2ckf2HYy8eSSjbh5F8f5i5vxqDkufWMrOLw9sKJzdLZtJqybVVZESxswWu3tepHkN4oxbRETiq6LuP/ft2kd6ZjqDLxiMmTF74my6HdONLkd2oXFWY8bNGEenwzsBkNYoje/87jt0GNShXlyPjgedcYugz5NITU3pPiVyI+Bu2bTu2ZqMZhmc8/w5AGz4YAPZOdk0zmp8wPLhEjGAR7Kot2fc7q4BIqTGUuHHq0gyW3TfoohJG0L3Tl/x4RVlGggfNOCgmLZbH65Hx0PKtipv0qQJmzdv1peu1EjJ7XNNmjRJdCgiKWPFMyv468C/sm/XPgDSMtJo1DTyeWB2TvYBd/VIzcTtjNvMHga+C2xw94HBtMHAfUATYD9wubu/U53td+nShXXr1mkcZamxJk2a0KVLl0SHIVLnKquKLqnVXPv2Wp678Dl+8OQPODj3YJq0akKr7q0o2FxAZvNMhlw8hIwmGbomXUfiWVU+DbgbeCxs2u+AW939X2Z2cvB6VHU2npGRQY8ePWoao4hIg1RRY7KuI7oy/aTpHPubY+k/rj9ZnbJoc0gbigtDvZH1OLYHPY4t+/0bj9tuJbK4JW53n2dm3ctPBkr6Bc0GvojX/kVEJLpoA27MuX4OV316FW37tqVJdugSUqturUobllVE16TrRl03TpsE/NvM/kDo+vrwaAua2URgIkBOTk7dRCci0kBU1JgsPSOds549K+J8Sby6bpz2/4Br3L0rcA3wULQF3X2qu+e5e1779u3rLEARkfpo4V0LeXr806WvG7eMfCtWvMeSlpqr68Q9ASj55DwJDK3j/YuINAir563m2QnPlo6SVbirkD3b9lC8P/T6lL+ekpCxpKXm6jpxfwGMDJ4fB3xawbIiIhKmomEut6/bzhu3vcGuDbsA2PHFDj6f8znb14aGNT76uqM554VzSGsU+tqvLwNuNERx6znNzGYQajHeDvgauBn4GLiT0LX1PYRuBztwnMtyIvWcJiLSkJRvBQ6Q3jid4yYfx/CfDufLd79kat5UznruLPqO7Uvx/mIs3dRJVYqqqOe0lO3yVESkIYnWpWjj7MZct+063J3dm3bTvH3zBEQnta2ixJ2yPaeJiDQERYWh4S3z10RuBb53+14AzExJu4FQ4hYRSVIv/+xlHhz2IBC9tbdagTc8StwiIkli7dtreerMp9i/dz8AHQ/vyCEnHUJRYRGjJ49WK3ABUnx0MBGRVLZ/z34+eeETug7vSlbHLPZs28Oa+WvYsnILBw04iNyzv2nhrS5FpYQap4mIxEG0ATyK9hWxb+c+mrZpyuZPN3N3n7sZc+cYhl09jOKiYswMS1NL8Iau3o7HLSKSjKIN4FFcXMxrv3qNPmP7cMpfT6Ft77ZcNP8iOg/tDEBauq5eSuWUuEVEalm0ATzm3jiXEdeNoN232pVO7zq8a12HJylOiVtEpBbt2banwgE8hl6hnp6lZlQvIyJSQ+5e2if4h7M+jLqcbt2S2qDELSJSA9vXb+fegffywcwPABh41kCO/fWxunVL4kaJW0Skila9sYqP//kxAFkds2jfvz1NWjcBILN5JsfccIwG8JC40e1gIiIxKCwoJKNp6Cz60WMfZe/2vUxcPDHBUUl9pdvBREQqEe2+a4D//vW/zL1xLpNWTyKzRSanPnwqLQ5ukeCIpaFSVbmINHgl913nr84HD913/eyFz7LgzgUAdMrrxKAJg0q7Im3do/UB17BF6ooSt4g0eJHuuy4uLGbebfMA6Dy0Myf+6USatW2WiPBEylDiFpEGL9qQmQVbC+o4EpHKKXGLSINUsLWAZTOWARoyU1KLEreINEiL7lvE0+c8zdbPt2rITEkpalUuIg3Czq938tKPX+LwiYfT47ge5F2WR++Te9O6Z2ta92wNaMhMSQ1K3CJSb3mxs/PrnWR1zKJJdhO+eu+r0uvZTVs3pWnrpqXL5o7PVaKWlKDELSL11t+/93d2rN/BpYsupVGTRlyx4gqNdS0pL26J28weBr4LbHD3gWHTrwKuAIqAF9z95/GKQUTqr0gdpvQ8oSdLpi3hqGuOIq1RGkMuGcK+XfvAAUNJW+qFeJ5xTwPuBh4rmWBmxwKnAYPcfa+ZHRTH/YtIPVXSYUrJvdf5q/OZPXE2h192OAv+tIDOR3Sm+6ju9D21b4IjFal9cUvc7j7PzLqXm/z/gDvcfW+wzIZ47V9E6q9IHaYU7i5kxVMruOKjK2jXt12CIhOJv7q+HawP8G0zW2hmb5jZEdEWNLOJZrbIzBZt3LixDkMUkWQXrcOU/LX5StpS79V14m4EtAGOBK4F/mFmES86uftUd89z97z27dvXZYwiksTe+PUbEOVStTpMkYagrhP3OuBpD3kHKAb081hEKvT10q/ZvWk3AH2+24f+Z/SnUdOyV/rUYYo0FHWduJ8FjgUwsz5AJrCpjmMQkRSy44sd3D/kft7+09sAdDysIz+Y+QNOfeBUsrtlg0F2t2zGTh2r+7ClQai0cZqZ/QB4yd13mNkNwBDgN+7+biXrzQBGAe3MbB1wM/Aw8LCZLQf2ARPc3WtYBhGpZzZ/spm1b61l8AWDyeqUxRl/P4Meo3uUWUYdpkhDFUur8hvd/UkzOxo4Hvg9cC8wrKKV3P3sKLPOrVqIItLQLJiygOUzltNvXD8aZzWm/xn9Ex2SSNKIpaq8KPh7CjDV3V8gVMUtIlIrdm3YxeyJs9m4InQHyahbRnHFR1fQOKtxgiMTST6xJO71ZnY/cCbwopk1jnE9EZEKlV4pM1jx9ArWL1wPQPODmtPi4BYJjEwkecVSVf5DYAzwB3ffZmYdCd3KJSISk0jdk25euZkNyzbww6d+SPP2zZm0ehKZzVWZJ1KZSs+c3X038Bywy8xygAzgo3gHJiL1Q0n3pPmr88G/6Z5088ebadyyMUWFoatxStoisYmlVflVhFqEf03ovmsIddl/aBzjEpF6Ilr3pGvfWsukVZMSFJVI6oqlqvzHQF933xzvYESkfnH36N2TRpkuIhWLpZHZWkD/YSJSZU+Pf5q0RpG/ZtQ9qUj1xHLG/Tnwupm9AOwtmejuf4pbVCKSsrav306LDi1IS08jd3wumS0yWTZ9WZnqcnVPKlJ9sZxxrwFeIXTvdlbYQ0SkjI0rNnJX77t476H3AOhzSh/GTh3L2Klj1T2pSC2xWHscNbMWAO6+M64RRZCXl+eLFi2q692KSAy82Nn6v6206dUGd+eNW99g8AWDadW9VaJDE0lZZrbY3fMizav0jNvMBprZe8AHwAdmttjMBtR2kCKSml688kUeHv4we3fsxcwYdcsoJW2ROIrlGvdU4CfuPhfAzEYBDwDD4xiXiCSxbau30bRNUxpnNWbIJUPoclQX3YctUkdiucbdvCRpA7j760DzuEUkIklt51c7uaffPfznjv8A0HFIRwadNwhLswRHJtIwxNSq3MxuBB4PXp9LqKW5iNRT5bsoPfbXx3LQgIPoOKQjLTq04MQ/n0jvk3snOkyRBqnSxmlm1hq4FTg6mPQmcIu7b41zbKXUOE2k7pR0URp++1ZaozTcnUmrJtGyS8sERifSMFTUOK3SM+4gQV9d61GJSFKK1EVp8f5imrZrSlZn3QkqkmhRE7eZTXH3SWY2m1Df5GW4+6lxjUxEEiJaV6QFmwsw03VskUSr6Iy75Jr2H+oiEBFJvC2fbSE7Jzs0klc56qJUJDlETdzuvjj4+0bdhSMiibLi6RX844x/cMwNx/D2H99WF6UiSaqiqvJlRKgiL+HuGtZTJMW5O3u27aFp66b0OrEXo24ZxVE/OYp2fduVaVU+evJodVEqkiSitio3s24Vrejuq+MSUQRqVS4SH8+c9wybP9nMRW9dRFp6LN06iEhdqFar8pomZjN7GPgusMHdB5ab91NC187bu/ummuxHRKqmqLCItEZpmBl9T+vLzq/qfPgBEamBqD+xzWyHmW2P8NhhZttj2PY0YEyE7XYFTiA06piI1KEdX+zg/sH3s+xvywDof0Z/hl45VGfbIimkojPuGt2w6e7zzKx7hFl/Bn4OPFeT7YtI7NwdM6NFhxYcNPAgmrdXr8UiqaqiM+6Wwd82kR7V2ZmZnQasd/f3Y1h2opktMrNFGzdurM7uRAT45PlPeHDog+zbuQ9LM86YeQa9TuiV6LBEpJoquo/7b4SuUS8m1Lo8vOcFB3pWZUdm1gz4FaFq8kq5+1RCI5ORl5cX26DhInKAJq2bkJaRxu7Nu8lsoRG8RFJdRVXl3w3+9qilffUCegDvB70vdQHeNbOh7v5VLe1DpEEqMyhI12w6De1E56GdGXHtCHJG5HDR/IvU65lIPRHL6GCY2emEBhlx4E13f7aqO3L3ZcBBYdtcBeSpVblIzZQfFCR/TT7b129n51c7Gf6z4ZiZkrZIPVJpU1Iz+ytwGbAMWA5cZmb3xLDeDOBtoK+ZrTOzi2sarIgcKNKgIF7kbF+7XQlbpB6K5Yz7OKCfBz21mNmjwAeVreTuZ1cyv3ssAYpIxaINChJtuoiktlhu3lwJ5IS97hpME5EE2rtjL8tmLIs6+IcGBRGpnyq6HWy2mf0TyAJWmNnrZjYXWBFME5EEWvDnBTxz3jMM+/EwMppllJmnQUFE6q+Kqso1nKdIktm3cx+7N++mVbdWDL92OIeMOYTOQzvT4qAWGhREpIGIOshIMtEgIyKh3s8eHv4wxUXFXLLwEjU8E6nHqjXIiIgkh/179pPeOB0zY+TNI8lolqGkLdKAaWQBkSS2fd127j30XpY+vhSAQ8YcQrdjKhxxV0TqOSVukSTWomMLOg/trBbiIlIqlg5YRpjZK2b2iZl9bmb/M7PP6yI4kYboy3e/ZPpJ09m7Yy9p6Wmc/sTpdB/VPdFhiUiSiOUa90PANYQGGymKbzgisn/PfjZ9tIlt/9vGwYcenOhwRCTJxJK48939X3GPRKQB+/LdL/lqyVccdtFhdB3elSs/ub7NQWYAABo/SURBVJL0jPREhyUiSSiWxD3XzH4PPA3sLZno7u/GLSqReq7MaF452WR1ymLnlzvJPSeXRk0aKWmLSFSxJO5hwd/w+8mcUB/mIlJFB4zmtTqfXRt2MWbKGBo10R2aIlKxSr8l3P3YughEpKGINJrX/oL9vHn7mxw+8fAERSUiqSJq4jazc939CTP7SaT57v6n+IUlUj/t3bFXo3mJSI1UdMbdPPirAUVEasEXi7/giROfoHn75uzasOuA+bpXW0RiETVxu/v9wd9b6y4ckfqrff/29D6pNwcPPpjXb3q9THW5RvMSkVip5zSROPpi0RfMOnsWRYVFZDTN4PuPf5/hPx3O2Kljye6WDQbZ3bIZO3WsRvMSkZioCatIHG1bvY0189eQvzqfNoe0KZ2eOz5XiVpEqkWJW6SW5a/NZ9NHm+j1nV70H9ef3if1JqNZRqLDEpF6otLEbWaNgXFA9/Dl3f22+IUlkrpevPxFvnzvS67+7GoaNW6kpC0itSqWM+7ngHxCfZXvrWRZkQZp7469mBmZLTI56e6TKC4splFjVWiJSO2L5Zuli7uPiXskIimqsKCQB/IeIOfbOZz64Km06tYq0SGJSD0WS6vyt8ysyq1ozOxhM9tgZsvDpv3ezD4ys6Vm9oyZ6RtOUl5G0wzy/l8egyYMSnQoItIAxJK4jwYWm9nHQcJdZmZLY1hvGlD+TP0VYKC7Hwp8AvyyStGKJInt67fz2PGP8dWSrwA4ctKRdPt2twRHJSINQSxV5SdVZ8PuPs/Mupeb9nLYywXAGdXZtkhdKj+S1+jJo+l1Yi+2r9tO/pp8OgzukOgQRaQBiWWQkdVmNgj4djDpTXd/vxb2fREwM9pMM5sITATIycmphd2JVF2kkbxmT5zN2KljuXz55aQ1Uh9GIlK3Kv3WMbMfA9OBg4LHE2Z2VU12ambXA/uD7Ubk7lPdPc/d89q3b1+T3YlUW6SRvAp3FzLn+jlK2iKSELFUlV8MDHP3XQBm9lvgbeCu6uzQzC4AvguMdnevzjZE6opG8hKRZBPLKYMBRWGvi4JpVWZmY4CfA6e6++7qbEOkLmU2z4w4XSN5iUiixJK4HwEWmtktZnYLoUZlD1W2kpnNIHRm3tfM1pnZxcDdhIYJfcXMlpjZfdUPXST+Rlw3gvTM9DLTNJKXiCSSxVJbbWZDCN0WBqHGae/FNapy8vLyfNGiRXW5S2mg3J03bn0DSzdG3jgSiNyqXAOEiEg8mdlid8+LNC/qNW4za+nu282sDbAqeJTMa+PuW2o7UJFEMzO2fr6VtPQ03B0z00heIpJUKmqc9jdCjcgWA+Gn5Ra87hnHuETq1MqXVtJ+QHuyu2Zz6kOnkp6RXvlKIiIJEDVxu/t3g7896i4ckbpXsKWAJ3/4JAPPGsjYqWOVtEUkqcVyH/ecWKaJpJqCrQUANG3TlPNePo8xd2osHRFJflETt5k1Ca5vtzOz1mbWJnh0BzrXVYAi8fDle1/yl55/4aNnPwKgy5FdyGiqcbNFJPlVdI37R8AkoBOh69wl925vJ3Rbl0jKOmjAQfQ7ox8HH3pwokMREamSqGfc7n5ncH37Z+7e0917BI9B7q7ELSln8yebeea8Z9i/Zz/pmemc+sCptO7ZOtFhiYhUSSyDjNxlZgOB/kCTsOmPxTMwkdq2ZeUWVr60kk0fb6LDII3oJSKpqdLEbWY3A6MIJe4XCQ3z+R9AiVuSUniHKS27tGTQBYM47rbj6H1yb67+/GoaZzVOdIgiItUWS5enZwCjga/c/UJgEKCOmiUplQzDmb86Hxy2r93Om79+k8VTFwMoaYtIyoslcRe4ezGw38xaAhuArvENS6R6Ig3DCfDm7W8mIBoRkdoXy7Cei8ysFfAAodblOwkNHiKSdDQMp4jUd7E0Trs8eHqfmb0EtHT3pfENS6R6snOyQ9XkEaaLiNQHFXXAMqT8A2gDNAqeiyQFd+fdh95l/TvrGT15NBnNynakomE4RaQ+qeiM+48VzHPguFqORaRaCncVMu/X8+gxugenPXQagIbhFJF6K6bxuBNN43FLJFs/30qr7q2wNCN/TT5ZnbNIS4+lvaWISHKr1njcYSufH2m6OmCRRNr08SbuP+x+jpt8HEddc5SuYYtIgxFLq/Ijwp43IXRP97uoAxZJoLZ92jLyppHknq0qcBFpWGJpVX5V+Ovg1rC/xy0ikSg2fbSJF698kdOfOJ0WHVpw9HVHJzokEZE6V50LgruAHrUdiEhlivcXs/WzrWxbvS3RoYiIJEws17hnE2pFDpAO9AP+Ec+gREoU7i5k5Usr6Xd6Pw4aeBBXfnIl6RnpiQ5LRCRhYrnG/Yew5/uB1e6+rrKVzOxh4LvABncfGExrA8wEugOrgB+6+9YqxiwNyPzfz2febfO48uMraXNIGyVtEWnwKq0qd/c3gI8JDSzShlDyjsU0YEy5adcBc9y9NzAneC1ygJL+xkf8fATnv3Y+bQ5pk+CIRESSQyxV5ZcANwGvAQbcZWa3ufvDFa3n7vPMrHu5yacRGiIU4FHgdeAXVYpY6p3wYTizc7JpP6A9BZsKuPDNC8lomkH3kd0THaKISNKIpar8WuAwd98MYGZtgbeAChN3FAe7+5fB86+Ag6MtaGYTgYkAOTk51diVpIKSYThLzrDzV+ez86ud9DyhJ5ZmCY5ORCT5xNKqfDOwI+z1jmBajXioy7ao3ba5+1R3z3P3vPbt29d0d5KkIg3DWbS3iA1LN5DWSL2giYiUF8sZ90pgoZk9RyjRngYsNbOfALj7n6qwv6/NrKO7f2lmHQmN7S0NmIbhFBGpmlhOaT4DnuWbs+PngP8BWcGjKv4JTAieTwi2JQ3UlpVbaNmlZcR56sJURCSyWHpOuxXAzFoEr3fGsmEzm0GoIVo7M1sH3AzcAfzDzC4GVgM/rF7Ykup2b97N1LypdBraiYLNBWWqyzUMp4hIdLG0Kh8IPE7oVjDMbBNwvrt/UNF67n52lFn6Rm7A3B0zo1nbZoyZMoYeo3uwZt4aDcMpIhKjSof1NLO3gOvdfW7wehRwu7sPj394IRrWs37Ytmobs86exSn3nkKHwR0SHY6ISNKqaFjPWK5xNy9J2gDu/jrQvJZikwYko3kGhbsL2b1pd6JDERFJWbEk7s/N7EYz6x48bgA+j3dgUj/s37OfRfcvwt1p3r45P3rvR/Q8vmeiwxIRSVmxJO6LgPbA08AsoF0wTaRSy/++nBcue4G189cCqFMVEZEaito4zcyaAJcBhwDLgJ+6e2G05UXCFWwpoGmbpgyaMIh2/drRZViXRIckIlIvVHTG/SiQRyhpnwT8vk4ikpQ37zfzuPfQe9m9eTdmpqQtIlKLKrodrL+75wKY2UPAO3UTkqS63qf0Zt+ufTRu2TjRoYiI1DsVJe7SanF332+ma5MS3fKZy8lfnc+In4+g42Ed6XhYx0SHJCJSL1VUVT7IzLYHjx3AoSXPzWx7XQUoqWHliyv5ZPYnFO8vTnQoIiL1WtQzbndPr8tAJDWEj52d1TGLEdeNYNhVwzjl3lNIy0jTiF4iInEWy+hgIsCBY2fv+GIH/77m3zRr00xdlIqI1BGdHknMIo2d7UXOnOvnJCgiEZGGR4lbYqaxs0VEEk+JW2LWsqvGzhYRSTQlbqnQnm17eOXnr7B/z36Ov/14MppllJmvsbNFROqWErdUaP0761kwZQGr31xN7vhcxk4dS3a3bDDI7pbN2Klj1TBNRKQOVToedzLQeNx1y4udr5d9TYdBoTGz89fmk91V1eEiInWlpuNxSwPz2o2v8dBRD5G/NtToTElbRCR56D5uKeXFjqUZQ68cStvebWnZJXJjNBERSRydcQsAc341h6fPfRp3J6tjFoMvGIz6pxcRST5K3AJAZotMGmc3xouSv82DiEhDlpCqcjO7BrgEcELjfV/o7nsSEUtDtvzvy2nVoxVdhnXh6F8erTNsEZEUUOdn3GbWGbgayHP3gUA6cFZdx9HQFe4u5NXrXmXhnQsBlLRFRFJEohqnNQKamlkh0Az4IkFxNDibP9lMm0PakNEsgwlzJ6gBmohIiqnzxO3u683sD8AaoAB42d1frus4GorwYThbdGjBrg27OOEPJ3DkpCNp3aN1osMTEZEqSkRVeWvgNKAH0AlobmbnRlhuopktMrNFGzdurOsw64WSYTjzV+eDw84vd2JpRnpjDbUuIpKqEtGq/Hjgf+6+0d0LgaeB4eUXcvep7p7n7nnt27ev8yDrg0jDcBYXFjP/t/MTFJGIiNRUIhL3GuBIM2tmoRZRo4EVCYij3tMwnCIi9U8irnEvNLOngHeB/cB7wNS6jqM+++r9r1g1dxXZOdmhavJyNAyniEjqSkgHLO5+s7t/y90Huvt57r43EXHUV+899B7zfzefb9/wbQ3DKSJSz6iv8noif00+xfuLad2zNcffcTwjbxpJs3bNyGyaWdqqPDsnm9GTR2sYThGRFKbEXQ8U7y9m2qhptOnVhvNeOY+MZhmlZ9q543OVqEVE6hEl7hS2b9c+MptnktYojbFTx9K6p+7LFhGp7zTISIra8tkW7vnWPSybsQyAnsf3VOIWEWkAlLhTVKturehxXA/a9mmb6FBERKQOKXGnkK+Xfs3M789k3659pDVK43uPfo9Oh3dKdFgiIlKHlLhTSMHWAta/s54tK7ckOhQREUkQJe4kl78mnw9nfQhA95Hdufqzq+kwqEOCoxIRkURRq/IkEz6aV3ZONlmds9jyyRYOOfEQMltk0qiJDpmISEOmLJBESkbzKhkYJH91Prs27GL05NFktshMcHQiIpIMVFWeRCKN5rW/YD8L7lyQoIhERCTZKHEnib079mo0LxERqZQSdxJY+/ZapuRMoflBzSPO12heIiJSQok7gbzYATj40IPpe2pfRvx8hEbzEhGRCqlxWoIsmLKAFbNWMOH1CWQ2z+R7j34PgBYHt9BoXiIiEpUSd4I0P7g5Lbu2pHB3IY2zGpdO12heIiJSEVWV15HCgkJeuPwFlk5fCkDu2bmM+9u4MklbRESkMjrjriPpmelsWLaBFh1aJDoUERFJYTrjjqPt67fz/GXPs2/nPtLS0zj/tfMZedPIRIclIiIpTIk7jrb9bxtLH1/K+nfWA5CekZ7giEREJNUpcdeyjR9uLL2OnXN0DpPWTKLHcT0SHJWIiNQXusZdQwcMCtIpi/zV+fQf159GTRrRrG2zRIcoIiL1SEISt5m1Ah4EBgIOXOTub8d7v+WTbE3vkY42KMgJfzxBo3iJiEhcJCq73Am85O5nmFkmEPfT0khJdvbE2QD0GN2D4qJiWnZuCcBHz35EZlYmPUf3BODla1+mTa825F2WB8CDwx6k26hufDDzg4iDgsz/7XyO+H9HxLtIIiLSANV54jazbOAY4AIAd98H7Iv3fiONvFW4u5A518+h6R+b0rJzS86efTYAc2+aS5tD2pQm7i/e+YLiwuLS9boM70K7b7XToCAiIlLnEnHG3QPYCDxiZoOAxcCP3X1X+EJmNhGYCJCTk1PjnVaUZE+66yQym38z3vU5L5xT5vUFb1xQZp0xfx4DwBu3vkH+6gO3q0FBREQkXhLRqrwRMAS4190PA3YB15VfyN2nunueu+e1b9++xjuNlkyzc7LpO7ZvmZbf2V2zadqmaaXbHD15tAYFERGROpWIxL0OWOfuC4PXTxFK5HEVjySbOz6XsVPHkt0tGwyyu2UzdupY9TUuIiJxU+dV5e7+lZmtNbO+7v4xMBr4MN77LUmmtT3ylgYFERGRupSoVuVXAdODFuWfAxfWxU6VZEVEJNUlJHG7+xIgLxH7FhERSWXq8lRERCSFKHGLiIikECVuERGRFKLELSIikkLM3RMdQ6XMbCOwuhY32Q7YVIvbSySVJfnUl3KAypKs6ktZ6ks5oPbL0s3dI/Y+lhKJu7aZ2SJ3rxet2lWW5FNfygEqS7KqL2WpL+WAui2LqspFRERSiBK3iIhICmmoiXtqogOoRSpL8qkv5QCVJVnVl7LUl3JAHZalQV7jFhERSVUN9YxbREQkJSlxi4iIpJCUTNxm1tXM5prZh2b2gZn9OJjexsxeMbNPg7+tg+lmZn8xs5VmttTMhoRta0Kw/KdmNiHK/iJuN1nKYWaDzeztYBtLzezMKPu7wMw2mtmS4HFJbZSjNssSzCsKi/GfUfbX2MxmBusvNLPuyVYWMzs2rBxLzGyPmX0vwv6S6bh8K/gs7TWzn5Xb1hgz+zgo53VR9heX41Jb5Yi2nQj7G2Vm+WHH5KbaKEdtliWYt8rMlgUxLoqyv6j/a8lSFjPrW+5/ZbuZTYqwv7gcl2qUY3zwXi4zs7fMbFDYtuL/f+LuKfcAOgJDgudZwCdAf+B3wHXB9OuA3wbPTwb+BRhwJLAwmN6G0LCibYDWwfPWEfYXcbtJVI4+QO/geSfgS6BVhP1dANydzMckmLczhv1dDtwXPD8LmJmMZQnbZhtgC9AsyY/LQcARwGTgZ2HbSQc+A3oCmcD7QP+6Oi61WI6I24mwv1HA88l8TIJ5q4B2leyv0s9nMpSl3GftK0IdkNTJcalGOYYT5ArgJL75Lq6T/5Na/1Am4gE8B3wH+BjoGHYgPg6e3w+cHbb8x8H8s4H7w6aXWa788uW3myzliLCd9wkSebnpFxCnBFGbZSG2xP1v4KjgeSNCPRZZspUlbNpEYHqU7SfNcQlb7hbKJryjgH+Hvf4l8MtEHZfqliPadiJMH0WcEndtloXYEndM3xuJLkvYvBOA+VHm1clxibUcwfTWwPrgeZ38n6RkVXm4oIrhMGAhcLC7fxnM+go4OHjeGVgbttq6YFq06eVF226tqWE5wrczlNAvvc+i7GpcUMXzlJl1rZ3oy6qFsjQxs0VmtsAiVC2XX9/d9wP5QNvaKkOJ2jouhH5Vz6hgV8lyXKKJ9X8l7selhuWItp1IjjKz983sX2Y2oLrxViGG6pTFgZfNbLGZTYyyTKzHrkZq67hQ+f9KXI9LNcpxMaEaDaij/5OUTtxm1gKYBUxy9+3h8zz0U6bW73WLx3Zrqxxm1hF4HLjQ3YsjLDIb6O7uhwKvAI/WKPDIMdRGWbp5qOvAc4ApZtartuOMRS0fl1xCv7IjSZXjknC1eEyibifwLqHP4SDgLuDZGgVexRiqUJaj3X0IoeraK8zsmNqOMxa1eFwygVOBJ6MsEtfjUtVymNmxhBL3L2ozjsqkbOI2swxCb/B0d386mPx18CVZ8mW5IZi+Hgg/i+kSTIs2vbxo202WcmBmLYEXgOvdfUGkfbn7ZnffG7x8EDi8tspRm2Vx95K/nwOvE/r1W17p+mbWCMgGNidbWQI/BJ5x98JI+0qy4xJNrP8rcTsutVSOaNspw923u/vO4PmLQIaZtauFYlQUQ5XLEva/sgF4BhgaYbFYj1211FZZAicB77r715FmxvO4VLUcZnYoof/X09y95DNeJ/8nKZm4zcyAh4AV7v6nsFn/BCYEzycQuk5RMv18CzkSyA+qP/4NnGBmrYPWgicQ+awo2naTohzBr9RngMfc/akK9tcx7OWpwIraKEew7doqS2szaxxssx0wAvgwwi7Dt3sG8FrwizhpyhK23tlUUPWXZMclmv8Cvc2sR/B5OyvYRnlxOS61VY4KtlN+uQ7BsiWXn9KovR8gtVWW5maWVfKc0PfX8giLVvb5rLZa/HyVqOx/JS7HparlMLMc4GngPHf/JGz5uvk/ifVieDI9gKMJVVksBZYEj5MJXSOYA3wKvAq0CZY34B5C132XAXlh27oIWBk8Lgyb/mDJctG2myzlAM4FCsO2sQQYHMy7DTg1eP5/wAeEGq/NBb6VbMeEUGvNZUGMy4CLw/YRXpYmhKrTVgLvAD2TrSzBvO6Efl2nldtHsh6XDoSuy20HtgXPWwbzTibU2vYzQjU7dXZcaqsc0bYTrHMZcFnw/MqwY7IAGJ5sx4RQy+X3g8cH5Y5JeFmifj6TpSzBvOaEknB2uX3E/bhUoxwPAlvDll0Utq24/5+oy1MREZEUkpJV5SIiIg2VEreIiEgKUeIWERFJIUrcIiIiKUSJW0REJIUocYukMDNrZWaXh73uZGZR7+Wv4b6+Z7UwGpOZ5ZrZtFoISaRB0u1gIinMQv0qP+/uA+tgX28Rug91U4zLN/JQP8yR5r0KXOTua2ozRpGGQGfcIqntDqCXhcYm/r2ZdTez5VA6zvezFhpHeJWZXWlmPzGz9yw0gEubYLleZvaShQaqeNPMvlV+J2bWB9jr7pvMLMvM/hd0EYmZtSx5bWavm9kUC40N/WMz+4GZLbfQoBDzwjY5m1CvUiJSRUrcIqntOuAzdx/s7tdGmD8QOJ1vxkDe7e6HAW8D5wfLTAWucvfDgZ8Bf42wnRGEBnjA3XcQ6kP+lGDeWcDT/k1f7JnunufufwRuAk700KAQp4ZtbxHw7WqUV6TBU+IWqd/muvsOd99IaOjA2cH0ZUB3C42GNBx40syWEBq7uWOE7XQENoa9fhC4MHh+IfBI2LyZYc/nA9PM7FIgPWz6BqBT9Yok0rA1SnQAIhJXe8OeF4e9Lib0/58GbHP3wZVsp4DQCEYAuPv8oFp+FJDu7uGDW+wKW+4yMxtG6Ox8sZkd7qGRlJoE2xSRKtIZt0hq2wFkVXdlD405/D8z+wGERkkys0ERFl0BHFJu2mPA3yh7tl2GmfVy94XufhOhM/aSIQ/7EHkkKxGphBK3SAoLzl7nBw3Afl/NzYwHLjazklGmTouwzDzgsJIhFQPTgdZUMAwj8HszWxY0mHuL0KhOAMcSGj9eRKpIt4OJSEzM7E5gtru/Grw+AzjN3c+r4nYaA28AR0e7XUxEotM1bhGJ1e3AMAAzuws4idDYw1WVA1ynpC1SPTrjFhERSSG6xi0iIpJClLhFRERSiBK3iIhIClHiFhERSSFK3CIiIink/wNlE3nQhALygwAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "fig = plt.figure(figsize=(8,4))\n", "plt.plot(t,w,'o:',color='purple',label='2 step Adams Method ')\n", "plt.title('Non Linear Population Equation')\n", "plt.legend(loc='best')\n", "plt.xlabel('time (yrs)')\n", "plt.ylabel('Population in billions')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": { "id": "EGPLOaZTaNb4" }, "source": [ "## Table\n", "The table below shows the time and the numerical approximation, $w$, for the non-linear population equation:" ] }, { "cell_type": "code", "execution_count": 24, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 696 }, "id": "RFRrAxT-aNb4", "outputId": "2ebd0e5c-2480-4923-8ffd-8b601c5a114c" }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
time t_iAdams approx of non-linear w
02000.06.000000
12001.06.084000
22002.06.933974
32003.07.869642
42004.08.848568
52005.09.851373
62006.010.857671
72007.011.846747
82008.012.799268
92009.013.698782
102010.014.532747
112011.015.292965
122012.015.975461
132013.016.579947
142014.017.109042
152015.017.567443
162016.017.961143
172017.018.296777
182018.018.581128
192019.018.820774
202020.019.021862
\n", "
" ], "text/plain": [ " time t_i Adams approx of non-linear w\n", "0 2000.0 6.000000\n", "1 2001.0 6.084000\n", "2 2002.0 6.933974\n", "3 2003.0 7.869642\n", "4 2004.0 8.848568\n", "5 2005.0 9.851373\n", "6 2006.0 10.857671\n", "7 2007.0 11.846747\n", "8 2008.0 12.799268\n", "9 2009.0 13.698782\n", "10 2010.0 14.532747\n", "11 2011.0 15.292965\n", "12 2012.0 15.975461\n", "13 2013.0 16.579947\n", "14 2014.0 17.109042\n", "15 2015.0 17.567443\n", "16 2016.0 17.961143\n", "17 2017.0 18.296777\n", "18 2018.0 18.581128\n", "19 2019.0 18.820774\n", "20 2020.0 19.021862" ] }, "execution_count": 24, "metadata": {}, "output_type": "execute_result" } ], "source": [ "\n", "d = {'time t_i': t, 'Adams approx of non-linear w': w}\n", "df = pd.DataFrame(data=d)\n", "df" ] }, { "cell_type": "markdown", "metadata": { "id": "224Fv07NaNb5" }, "source": [ "## 3. Non-Linear Population Equation with an oscilation \n", "\\begin{equation} y^{'}=0.2y-0.01y^2+\\sin(2\\pi t), \\ \\ (2000 \\leq t \\leq 2020), \\end{equation}\n", "with the initial condition,\n", "\\begin{equation}y(2000)=6.\\end{equation}\n", "\n", "## Specific 2 Step Adams Bashforth for the Non-Linear Population Equation with an oscilation\n", "To write the specific \n", "\n", "\\begin{equation} w_{i+1}=w_{i}+\\frac{0.1}{2} \\big[ 3(0.2w_{i}-0.01w_{i}^2+\\sin(2\\pi t_{i}))- (0.2w_{i-1}-0.01w_{i-1}^2+\\sin(2\\pi t_{i-1}))\\big] \\end{equation}\n", " \n", "for $i=1,...,199$, where $w_i$ is the numerical approximation of $y$ at time $t_i$, with step size $h$ and the initial condition\n", "\\begin{equation}w_0=6.\\end{equation}\n", " As $w_1$ is required for the method but unknown we will use the numerical solution of a one step method to approximate the value. Here, we use the 2nd order Runge Kutta approximation (see [Runge Kutta notebook](https://github.com/john-s-butler-dit/Numerical-Analysis-Python/blob/master/Chapter%2003%20-%20Runge%20Kutta/01_2nd%20Order%20Runge%20Kutta%20Population%20Equations.ipynb) )\n", "\\begin{equation}w_1=6.11.\\end{equation}" ] }, { "cell_type": "code", "execution_count": 25, "metadata": { "id": "ZwZH25t9aNb6" }, "outputs": [], "source": [ "def nonlin_oscfun(t,w):\n", " ftw=0.2*w-0.01*w*w+np.sin(2*np.math.pi*t)\n", " return ftw" ] }, { "cell_type": "code", "execution_count": 26, "metadata": { "id": "aAPnbMp0aNb7" }, "outputs": [], "source": [ "## INSERT METHOD HERE\n", "w=np.zeros(N+1)\n", "w[0]=6\n", "w[1]=6.11\n", "for n in range(1,N):\n", " w[n+1]=w[n]+h/2*(3*nonlin_oscfun(t[n],w[n])\n", " -nonlin_oscfun(t[n-1],w[n-1]))" ] }, { "cell_type": "markdown", "metadata": { "id": "tccu5BOqaNb8" }, "source": [ "## Results\n", "The plot below shows the numerical approximation, $w$ (circles) for the non-linear population equation:" ] }, { "cell_type": "code", "execution_count": 27, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 295 }, "id": "rcg4LHuxaNb8", "outputId": "0208e444-1185-47f9-a8ac-309f843d9064" }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAe4AAAEWCAYAAACg1nQiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXxU5dn/8c+VEAJhCTvIEoIoguwY3HBBaakbitVWLaK48VSrFq22+vj83J7qYyutVq0LdUGrUq077ooirihQFgU3NkVF9gCyJrl+f5yTOISZZLLMTCb5vl+veWXmLPe57jmTuebc55z7NndHRERE0kNGqgMQERGR+Clxi4iIpBElbhERkTSixC0iIpJGlLhFRETSiBK3iIhIGlHilpQys8lm9scarL/ZzPaszZjqGjP7xMyG19Xtm9l0Mzs3iSGlnJkNN7MVSdxe2T4ws2vN7OHaKEvSkxK3lDGzZWa2NUyG34dJtXmq4yoVLUG4e3N3X5KAbUW+F6WPO2p7O1G2u9sPGXfv6+7TE73tWCK3X9OkIdVT3c9AXfw8Sc0pcUt5o9y9OTAEKAD+J8XxpNKo8IdB6ePCVAckIqLELVG5+zfAS0A/ADM7Pmxi2xAe+fYpXTY8Or3SzBaa2Xoze8DMmoTzxpnZO5Flm5mb2V7lt2lmrc3seTNbHZbzvJl1DefdABwK3BF59BtZlpnlmtlD4frLzex/zCwjMg4zmxiWvdTMjq7Oe2NmmWE5a8xsiZn9JoyjUcT78ZOI5Xc5SjWzf5vZSjMrNLMZZtY3nD4eGAP8Pqzj1PLlmVm2md1qZt+Gj1vNLDucN9zMVpjZ78xslZl9Z2ZnxajDEWa2IOL1a2b2UcTrt81sdOT2zewo4L+BU8L45kUU2d3M3jWzTWb2qpm1i7HdduF+3WBm68LtlO6jzmb2ZLj/lprZxRHr7W9m74frfWdmd5hZ43CemdktYZ03mtkCMyv93Fb7M2FmZ5nZorBOS8zsv2J+KHav58Fm9lG4jz8ys4Mj5o0Ly9sUbnNMxLzzIra50MyGRO6DGNtK+edJkkuJW6Iys27AMcB/zKwXMAWYALQHXgSmln5xhsYAPwN6Ar2o3pF6BvAA0B3IA7YCdwC4+1XA28CFFRz93g7kAnsChwNnAJFfNAcAnwHtgD8D95mZVSPO84DjgMEErRInV3H9l4C9gQ7AHOARAHefFD7/c1jHUVHWvQo4EBgEDAT2Z9f3uhPBe9AFOAf4u5m1jlLOB8DeYSLNAgYAnc2shZk1Dev1duQK7v4ycCPwWBjfwIjZvyJ4rzsAjYHLYtT9d8AKgs9RR4IfAh4m06nAvDD2EcAEM/tZuF4xcAnBvjsonH9BOG8kcBjB5y4X+CWwNpxXk8/EKoL93DJc55bSRFoRM2sDvADcBrQF/gq8YGZtzaxZOP1od28BHAzMDdf7BXBtGGNL4PiIelSkLnyeJImUuKW8Z8xsA/AO8BbBF/UpwAvu/pq77wQmAk0JvnRK3eHuX7v7OuAG4LSqbtjd17r7k+6+xd03heUcHs+6ZpYJnApc6e6b3H0Z8BdgbMRiy939H+5eDDwI7EGQPGJ5JjzCK32cF07/JXBrRH3/r4r1vD+McTvBF/VAM8uNc/UxwPXuvsrdVwPXsWsdd4bzd7r7i8BmYJ8oMWwFPiJIePsRJMx3gWEEX+RfuHs8SaPUA+7+eVju4wSJIJqdBO979zDGtz0YMGEo0N7dr3f3HeF1C/8g2Ke4+2x3/8Ddi8J9ew8/fjZ2Ai2A3oC5+yJ3/66mnwl3f8HdF3vgLeBVglafyhxL8P79M4x3CvApUJo4S4B+ZtbU3b9z90/C6ecSJNmPwm1+6e7LK9tYXfg8SXIpcUt5o929lbt3d/cLwi/izkDZF4i7lwBfE/wKL/V1xPPl4TpVYmY5ZnZP2KS5EZgBtAq/gCvTDsiKjDN8HhnjytIn7r4lfFrRxXel70Xp4x/h9M7sXt+4WNDMfpOZLQ7ruCwi/njssi/Y/b1e6+5FEa+3ELuObwHDCZL3W8B0gmR4ePi6KlZGPK9omzcDXwKvhs3FV4TTuxMc8Zf9UCI4Gu8IYGa9LGhiXxm+bzcSvmfu/gZBy8zfgVVmNsnMWlLDz4SZHW1mH1jQpL+BoAUqnv1Ufh+VbdfdfyD4Ifxr4Dsze8HMeofLdAMWx1F+mTr2eZIkUeKWeHxL8MUKBOcUCb5kvolYplvE87xwHYAfgJyIdTtVsJ3fEfyaP8DdWxIkFIDSpsuKhrJbQ3B00D1iWl65GGvLd+xe30i71JmgubHUr4ATgJ8QNEHmh9PjqSOU2xfs+l5XVfnE/RaVJ+4aDScYHhn+zt33JGgKvtTMRhD8EFpa7odSC3c/Jlz1LoKj1r3Dz8Z/8+N7hrvf5u77AfsSNJlfTg0+E+F53icJWpc6unsrglNE8ZxaKb+Pdtmuu7/i7j8lOLr/lKBlgfA96BlH+ZHq0udJkkSJW+LxOHCsmY0Iz4f+DtgOvBexzG/MrGt4fu8q4LFw+jygr5kNsuCCtWsr2E4LgvPaG8Jyrik3/3uCc5W7CZs6HwduCM/TdgcuBRJx69LjwMVhfVsDV5SbPxc41cyyzKz8OfAWBO/dWoLkfmO5dWPWMTQF+B8za2/BBWBXU/06vkfwQ2l/4MOwybY7wXnfGTHW+R7ID89JV5mZHWdme4U//goJzl2XAB8Cm8zsD2bWNDyS7GdmQ8NVWwAbgc3hEer5EWUONbMDws/mD8A2oKSGn4nGQDawGiiy4KK1kXFW80Wgl5n9yswamdkpBD8onjezjmZ2QniueztB03NJuN69wGVmtp8F9gpjrkhd+jxJkihxS6Xc/TPgdIILfdYQnKsb5e47IhZ7lOAc4BKC5r4/hut+DlwPvA58QXDuPJZbCc6dryG4eOrlcvP/BpxswRXAt0VZ/yKCL+4l4XYeBe6Pu6K7m2q73sf9dDj9H8ArBD9K5gBPlVvv/xEcOa0nOGf4aMS8hwiaI78BFhLUM9J9wL5hc/EzUWL6IzALmA8sCLdfrQ5swmbbOcAnEfvyfYLzvqtirPbv8O9aM5tTjc3uTfBZ2Bxu6053fzNMsscRnBtfSvAZuJfgKBKCi91+BWwieP8fiyizZThtPcF7u5agSR6q+ZkIr7G4mCDxrw+3/Vw8FQyvDTiO4AfuWuD3wHHuvobgO/dSgqPadQStG+eH6/2b4LqOR8N6PgO0qWRzdebzJMljwXUhItVnZsuAc9399VTHkgpmlk+QbLLKnQ8UEal1OuIWERFJI0rcIiIiaURN5SIiImlER9wiIiJppFGqA4hHu3btPD8/P9VhiIiIJMXs2bPXuHv7aPPSInHn5+cza9asVIchIiKSFGYWs0dGNZWLiIikESVuERGRNKLELSIikkbS4hx3NDt37mTFihVs27Yt1aFIgjRp0oSuXbuSlZWV6lBEROqMtE3cK1asoEWLFuTn5/PjuPdSX7g7a9euZcWKFfTo0SPV4YiI1Blpm7i3bdumpF2PmRlt27Zl9erVqQ5FRCSmBY8sYNpV0yj8qpDcvFxG3DCC/mP6J3SbaZu4ASXtek77V0TqsgWPLGDq+Kns3LITgMLlhUwdPxUgoclbF6eJiIhUw7SrppUl7VI7t+xk2lXTErpdJe4aeuaZZzAzPv3006jzhw8fnrTOY8aNG0ePHj0YNGgQvXv35rrrrqtWOddeey0TJ06MOu+2226jT58+jBkzJu7y5s6dy4svvhhX+SIiibTgkQXcmn8r12Vcx635t7LgkQUVLv/Dqh9Yt3hd2etnz36W6ddOB6Dwq8Ko68SaXlsaTOKu6s6K15QpUzjkkEOYMmVKrZRXUzfffDNz585l7ty5PPjggyxdurRWy7/zzjt57bXXeOSRR+JavqioaLfELSKSCqVN24XLC8F/bNqe9+C8smXm3DeHd29+t+z1lOOn8Pz458telxSVUFJcAkBuXm7U7cSaXlsaROKOtbNqmrw3b97MO++8w3333ce//vUvALZu3cqpp55Knz59OPHEE9m6dWvZ8ueffz4FBQX07duXa665pmx6fn4+V155JYMGDaKgoIA5c+bws5/9jJ49e3L33XcD8N1333HYYYcxaNAg+vXrx9tvv11hbKW3yTVr1gyA66+/nqFDh9KvXz/Gjx9P6ahwt912G/vuuy8DBgzg1FNPLVt/4cKFDB8+nD333JPbbrsNgF//+tcsWbKEo48+mltuuYV169YxevRoBgwYwIEHHsj8+fOB4Ih67NixDBs2jLFjx3L11Vfz2GOPMWjQIB577LGY5YuIJFKspu3nzn2u7PXSaUv5fOrnZa+PuP4IDrv6sLLXJz50Ikf+75EAjLhhBFk5u96umpWTxYgbRiQi/DL1JnFPHj6ZuZPnAlC8s5jJwycz/+Egkbx+5etRd9bLE14GYMuaLUwePpnPpn4GwOaVm+Pa5rPPPstRRx1Fr169aNu2LbNnz+auu+4iJyeHRYsWcd111zF79uyy5W+44QZmzZrF/Pnzeeutt8oSHUBeXh5z587l0EMPZdy4cTzxxBN88MEHZQn+0Ucf5Wc/+xlz585l3rx5DBo0KGpMl19+OYMGDaJr166ceuqpdOjQAYALL7yQjz76iI8//pitW7fy/PPBL8ibbrqJ//znP8yfP7/sRwLAp59+yiuvvMKHH37Iddddx86dO7n77rvp3Lkzb775JpdccgnXXHMNgwcPZv78+dx4442cccYZZesvXLiQ119/nSlTpnD99ddzyimnMHfuXE455ZSY5YuI1Iai7UVlzxc9vYgpx0/B3WM2YZcUlZQ9//nDP+esGWeVve45sif5h+dHXa//mP6MmjSK3O65YJDbPZdRk0Yl/KryepO4K7Jxxcao07es3VKjcqdMmVJ2lHrqqacyZcoUZsyYwemnnw7AgAEDGDBgQNnyjz/+OEOGDGHw4MF88sknLFy4sGze8ccfD0D//v054IADaNGiBe3btyc7O5sNGzYwdOhQHnjgAa699loWLFhAixYtosZU2lS+cuVKpk2bxnvvvQfAm2++yQEHHED//v154403+OSTT8piHDNmDA8//DCNGv14k8Gxxx5LdnY27dq1o0OHDnz//fe7beudd95h7NixABx55JGsXbuWjRs3ltWnadOmMd+7eMoXkYYtnlOcW9dt5bOpn7Hjhx0AzJ08lxub3cjm74MDsO0bt7Ppm01sW78tdtN29x+nW0bV7mbpP6Y/E5ZN4JqSa5iwbELCkzak+e1gkcZNH1f2PDMrc5fXuXm5QTN5OaU7Maddzi7LN+/UvNLtrVu3jjfeeIMFCxZgZhQXF2NmDB48OOryS5cuZeLEiXz00Ue0bt2acePG7dLrW3Z2NgAZGRllz0tfFxUVcdhhhzFjxgxeeOEFxo0bx6WXXkqLFi3KLkC79957d9le8+bNGT58OO+88w5DhgzhggsuYNasWXTr1o1rr722bNsvvPACM2bMYOrUqdxwww0sWLBgl3gAMjMzKSoqoipKm+hjqWn5IlK/xbrVavOqzaz7Yh0H/PYA2u3Tjq/f+5p/Hf8vznrnLPKG5dFpcCcOufKQsttJB505iEFnBi2UI24YsUuZkJym7drWII64E3Ee4oknnmDs2LEsX76cZcuW8fXXX9OjRw/2228/Hn30UQA+/vjjsubwjRs30qxZM3Jzc/n+++956aWXqrS95cuX07FjR8477zzOPfdc5syZw4knnlh2IVpBQcEuyxcVFTFz5kx69uxZlqTbtWvH5s2beeKJJwAoKSnh66+/5ogjjuBPf/oThYWFbN4c32kCgEMPPbTsIrXp06fTrl07WrZsudtyLVq0YNOmTVWqr4g0XDs27+D1K6Kf4nxv4nsseHQB65esByDv0DzOfvds9hiyBwCdBnbiyP89kmYddj94SFXTdm2rN0fcFSndKbXZu82UKVP4wx/+sMu0k046if/85z9s3bqVPn360KdPH/bbbz8ABg4cyODBg+nduzfdunVj2LBhVdre9OnTufnmm8nKyqJ58+Y89NBDUZe7/PLL+eMf/8iOHTsYMWIEP//5zzEzzjvvPPr160enTp0YOnQoAMXFxZx++ukUFhbi7lx88cW0atUq7piuvfZazj77bAYMGEBOTg4PPvhg1OWOOOIIbrrpJgYNGsSVV15ZpXqLSP3m7nw/73syG2fSft/2bN+4nT+1+RNe7FGX3/zdZq4uvrrsiLpJbhO6Hdwt7u31H9M/7RJ1eVZ6dXFdVlBQ4OXvhV60aBF9+vRJUUSSLNrPIukj3u4/Fz29CDOj9+jeuDsTO06k17G9OOGBEwD44NYPePdP70a9UDi3ey4Tlk1IeF1Szcxmu3tBtHkN4ohbREQSq6LuP7dv3s7m7zYz/NrhALz/l/fJaJRB79G9MTN+8fgvaJX/Y2vfgRMOpFn7ZvXifHQiKHGLiEiNTfvv2N1/7vmTPVn35Y+9j/3i8V+Q0z6n7HX+8PzdykvEKc76Iq0Tt7trIIp6LB1O44g0VF7iYMFgQHMfnFth95+j/jFql+/qFp2j385aXn04H50IaXtVeZMmTVi7dq2+3Oup0vG4mzRpkupQRITgf7K0q8/Fry1mYseJrPsiOIrO7ZZLVrOsqOvl5uXqAKuWJeyI28zuB44DVrl7v3DaIOBuoAlQBFzg7h9Wp/yuXbuyYsUKjddcjzVp0oSuXbumOgyRequii8m8xCneWUyj7EasXrSayYdN5vj7j2efUfvQes/W7H3M3mUHTj2O7MGoe0bpnHSSJLKpfDJwBxB539Kfgevc/SUzOyZ8Pbw6hWdlZdGjR4+axigi0iBFu5jsufOCPrt7jerFbT1vY9gfhnHwZQfTukdreo3qRfOOQedUbXq2YfSDo3cpT+ekkydhidvdZ5hZfvnJQGkPHbnAt4navoiIxBZtwI2irUVMu2oa/cf0Z9BZg+g0qBMAjZo04oT7T6i0TJ2TTo5kX5w2AXjFzCYSnF8/ONaCZjYeGA/BABwiIlIzkRf0RusGGn4cS/qnf/5p0uKSqkn2xWnnA5e4ezfgEuC+WAu6+yR3L3D3gvbt2yctQBGR+mjmbTO5s++dwdXgQNO20QcBSvRY0lJzyU7cZwJPhc//Deyf5O2LiDQIy95axn0H31c2CmJuXi5dD+rKjs3BKFpH/+3olIwlLTWX7MT9LXB4+PxI4Iskb19EJG1VNMxl4deFPHfuc6ycuxKArKZZeLGz+bug29Deo3tzwn0nkN0yGJmvvgy40RAlrK9yM5tCcMV4O+B74BrgM+BvBOfWtxHcDja7srKi9VUuItKQlL8KHCAjK4ODLjuIn9z4E7as3cIdve7gmDuPod8p/VIYqdSGivoqT9tBRkREGpJb82+NekFZ4xaNuXJjMOpeSXEJGZlp26+WRKgocWsPi4jUUe5edpV3rC5FS89ZA0raDYT2sohIHfXa5a9x14C7KNpeFPNqb10F3vAocYuI1BHL317Onf3upPDr4Oi67y/7MvIvI8FhxA0jdBW4AErcIiIps61wG2//39t8OyvoRLJ5p+bktM1h69qtAHTZvwtDzhlCoyaNdBW4lEnrYT1FROqqWAN4FH5dyPbC7XTo14GMzAze/uPbWIbRuaAzbfduy7i3xsUsU12KCihxi4jUumgDeEwdPxXHeeu6t2iV34qxr46lcfPGTPhqAjltc1IcsaQTJW4RkVoWbQCPnVt28sZVbzB68mhadm1ZNl1JW6pKiVtEpBat+3JdhQN45A/PT25AUu/o4jQRkRravnE72zZsA2Dzys0xl9OtW1IblLhFRGpgW+E2/tr1r7x/y/sAdBvWjePuOU63bknCKHGLiFTRnHvnMP266QA0yW3C8OuG0/uE3gCYGfuN30+3bknC6By3iEgl3J3v5nxH5/06A/DNR9+w9rO1+NWOmXHQJQftto5u3ZJEUeIWESH2fdcAH97xIS9f/DIXfnYhbXu15ZjbjyGzcWaKI5aGSolbRBq8aPddP33m06xauIoRN4yg7y/60qRVE1p2C27jUtKWVNI5bhFp8KLdd+3Fzux7ZgNBV6QDxw4kq2lWtNVFkkqJW0QavFhDZm5dtzXJkYhUTolbRBqklXNXMvW/plJSVKIhMyWtKHGLSINRUlxC8Y5iADYs28DCfy9kzWdrNGSmpBUlbhFpELas2cLte9/OrHtmAdBrVC8u+foSOvTtoCEzJa3oqnIRqbc2fbuJVR+voufInuS0y2Gvo/ei3T7tAMjIzKBxs8Zly+q+a0kXStwiUm+9dPFLLJ+xnEu/uZTMrEyO/fuxqQ5JpMYSlrjN7H7gOGCVu/eLmH4R8BugGHjB3X+fqBhEpP6K1mFKhwEdmHbFNI6//3iad2zOiBtHkNEog8ws3Xct9Uciz3FPBo6KnGBmRwAnAAPdvS8wMYHbF5F6qrTDlMLlheBBhylTx09l8auLWTlvJWs/WwtA215tab1n6xRHK1K7EnbE7e4zzCy/3OTzgZvcfXu4zKpEbV9E6q9oHabs3LKTD2//kAnLJ5CRqetupf5K9qe7F3Comc00s7fMbGisBc1svJnNMrNZq1evTmKIIlKX7di8I2aHKYVfFSppS72X7E94I6ANcCBwOfC4mVm0Bd19krsXuHtB+/btkxmjiNRh7/3lPfDo89RhijQEyU7cK4CnPPAhUAK0S3IMIpJGirYXMfO2mXw761sADrjoAA6/5nB1mCINVqWJ28yamVlG+LyXmR1vZtXtaf8Z4IjSsoDGwJpqliUiDUDxjmJm/O8MFj65EICmbZoy/Nrh6jBFGixzj9HmVLqA2WzgUKA18C7wEbDD3cdUst4UYDjBEfX3wDXAP4H7gUHADuAyd3+jsiALCgp81qxZlS0mIvXEx499zKdPf8pJU07CzNj07SZadG6R6rBEksbMZrt7QbR58VxVbu6+xczOAe509z+b2dzKVnL302LMOj2ObYpIA1NSVIJlGJZhbFu/jQ3LNrBtwzaatm6qpC0SIZ5z3GZmBwFjgBfCaerNQERqzYZlG/h7n7+z6OlFAAw5bwjnvH8OTVs3TXFkInVPPIn7t8CVwNPu/omZ7Qm8mdiwRKS+KykuYcOyDQC07NaSToM70bRNkKgzMjOIccOJSINX6TnuukDnuEXSW7TuSb948Qu+fv9rLvzsQnVJKlJOjc5xh1d/XwbkRy7v7kfWVoAiUn+Vdk9a2tNZafekB15yIL2O76UOU0SqKJ6L0/4N3A3cSzAwiIhI3GJ1Tzr/4flMWDYhRVGJpK94EneRu9+V8EhEpN7ZuGJjhd2TikjVxdNGNdXMLjCzPcysTekj4ZGJSNp7ecLLMS8yU/ekItUTT+I+k6Bf8feA2eFDV4qJyG68xJn30Dx+WPUDACMnjmTkX0aqe1KRWlRp4nb3HlEeeyYjOBFJL+uXrue5c55jzr1zAGiV34oDJxyo7klFalE8XZ5mEYyjfVg4aTpwj7vvjLlSLdPtYCJ11/ol61n65lKGnDMEgG9nfcseQ/bAMnQftkh1VXQ7WDxN5XcB+wF3ho/9wmkiIsy8bSavXvoqW9dvBaBzQWclbZEEiueIe567D6xsWiLpiFuk7igpKmHOfXPIOySPDn07sHX9VnZu2UnLLi1THZpIvVHTI+5iM+sZUdie6H5ukQZrW+E2pl0xjXkPzQOgaeumStoiSRTPfdyXA2+a2RLAgO7AWQmNSkRSqnwXpfv/Zn92/LCDw685nJy2OZw36zxa79k61WGKNEiVJm53n2ZmewP7hJM+c/ftiQ1LRFIlWhel066ahmUag84aRKvurWjTU105iKRKzMRtZke6+xtm9vNys/YyM9z9qQTHJiIpEK2L0pKdJbTs2JJW3VulKCoRKVXREffhwBvAqCjzHFDiFqmHYnVFuvGbjUmORESiiZm43f2a8K/OZ4s0ANs3bWfWXbPI7ZYbNXmri1KRuqGipvJLK1rR3f9a++GISKosfWMpr//hdYb9YRgf3v7hLs3l6qJUpO6oqKm8RdKiEJGUKPyqkDWfraHnT3uyz/H7cMEnF9B+3/Z07N9xl6vKR9wwQl2UitQRlXbAUheoAxaRxHjkmEdYtWAVFy+5mMyszFSHIyKhijpgqaip/LaKCnX3iyvZ6P3AccAqd+9Xbt7vgIlAe3dfU1E5IlK7vvnwG9ru05YmuU04+vajyWiUoaQtkkYqaiqfXcOyJwN3AA9FTjSzbsBI4Ksali8iVbRxxUbuH3Y/B112ED/5v5/ofmyRNFTRVeUP1qRgd59hZvlRZt0C/B54tibli0h83J3VC1fToW8HWnZtyUn/OomeP+1Z+YoiUifF7KvczG4N/041s+fKP6qzMTM7AfjG3edVM14RqaL3//o+9wy+h7WfrwVg35P2JbtldoqjEpHqqqip/J/h34m1sSEzywH+m6CZPJ7lxwPjAfLy8mojBJEGo3hnMds3bienbQ4DzxhIo+xG6ltcpJ5I6FXlYVP58+7ez8z6A9OALeHsrsC3wP7uvrKicnRVuUjFdhkUpFsuGVkZtOvdjtOmnoaZxsYWSTfVGtbTzPY2s8lm9lcz62pmL5nZZjObZ2ZRC6uIuy9w9w7unu/u+cAKYEhlSVtEKlY6KEjh8kLw4N7sjV9vpO0+bZW0ReqhisbjfgB4j+CoeCZwP9AOuAz4e2UFm9kU4H1gHzNbYWbn1DxcESkv2qAgxTuKWfTkohRFJCKJVNE57ubuPgnAzH7t7v8Op79mZjdXVrC7n1bJ/Py4oxSRmGINChJruoikt4qOuEsinpcfFqgEEUmpL1/+kkePfZTcbtEH/9CgICL1U0WJu7eZzTezBRHPS1/vk6T4RCSGnVt3snHFRg667CCycrJ2madBQUTqr4qayvskLQoRicsXL37Bjh920PcXfelzYh/2GbUPGY0yyGmTo0FBRBqIinpOW57MQESkYu7Ou396Fy9x9j15X8yMjEZBo1n/Mf2VqEUaiIqOuEWkDljy+hK67N+F7JbZnPz4yTRp1US3eYk0YBWd4xaRFFu/dD0PH/Uw7018D4DmHZvTKFu/t0UaMn0DiNRBpeeqW/dozWlTT6PHET1SHZKI1BGVHnGb2VP/CLgAABpZSURBVDAze83MPjezJWa21MyWJCM4kYZo7uS53L737Xy/4HsA9j56bxo10W9sEQnE821wH3AJwfjcxYkNR6ThKikqIaNRBr1G9eLgxQfTZi+NlS0iu4sncRe6+0sJj0SkAdllUJC8XDr07wAOp009jZy2ORz5v0emOkQRqaPiSdxvhl2cPgVsL53o7nMSFpVIPVY6KEhp/+KFywvZ9O0m8o/Mp6SohMyszBRHKCJ1WTyJ+4Dwb+SIYA7okECkGqINClKys4S1n65V0haRSlWauN39iGQEItJQaFAQEamJmInbzE5394fN7NJo8939r4kLS6T+8RJn/iPzye2WGzVJa1AQEYlHRbeDNQv/tojxEJEqWPbWMp454xn2OmYvDQoiItVm7p7qGCpVUFDgs2bNSnUYItWyeeVmmndqDsDSN5eSPzyfjx/9WIOCiEhMZjbb3QuizlPiFkmc9295n7eue4vz55+vpnARiVtFiVvdMYkkUO/Rvdm6dmvZEbeISE1pkBGRWvbexPd44TcvANC6R2uO/OORZDbWbV4iUjsqPeI2s2zgJCA/cnl3vz5xYYmkry1rt7Bl1ZayLkxFRGpTPE3lzwKFBH2Vb69kWZEG6ZN/f0K7fdrRcUBHjvzjkViGacxsEUmIeBJ3V3c/KuGRiKSp7Zu28/LFL9NzZE9GPziajEwdZYtI4sTzDfOemVX5PhUzu9/MVpnZxxHTbjazT81svpk9bWatqlquSF2xetFq3J3sFtmc+eaZjLp3VKpDEpEGIJ7EfQgw28w+CxPuAjObH8d6k4HyR+qvAf3cfQDwOXBllaIVqSO++fAb7up/F/MenAdAu97t1M+4iCRFPE3lR1enYHefYWb55aa9GvHyA+Dk6pQtkky7DMHZLZcRN46g32n9GHHjCHqP7p3q8ESkgYlnkJHlZjYQODSc9La7z6uFbZ8NPBZrppmNB8YD5OXl1cLmRKputyE4vyrkufOeA2DY74elMjQRaaAqbSo3s98CjwAdwsfDZnZRTTZqZlcBRWG5Ubn7JHcvcPeC9u3b12RzItUWbQjOoq1FTLtqWooiEpGGLp6m8nOAA9z9BwAz+xPwPnB7dTZoZuOA44ARng79rUqDpiE4RaSuiefiNAOKI14Xh9OqzMyOAn4PHO/uW6pThkgylR/Fq5T6HReRVIkncT8AzDSza83sWoKLyu6rbCUzm0JwZL6Pma0ws3OAOwiGBH3NzOaa2d3VD10kMbZt2MbmlZsBOOrWo2jUdNeGKQ3BKSKpFM/FaX81s+kEt4UBnOXu/4ljvdOiTK404YukUklxCQ8c+gDN92jO2FfHMuTcIWQ1zdIQnCJSZ8Qc1tPMWrr7RjNrE22+u69LaGQRNKynJJq7l3VRuvDJhbTs2pKuB3RNcVQi0lBVd1jPRwkuIpsNRGZ3C1/vWWsRiqTQ1vVbeWrMUwy9YCi9juvFviftm+qQRERiipm43f248G+P5IUjknxZTbPYtmEbW9boekkRqfviuY97txtWo00TSSclxSV8dOdHFG0volGTRpz9ztkMGjco1WGJiFQq5hG3mTUBcoB2ZtaaH28Bawl0SUJsIgnz1Ttf8eJvXiQ7N5sBYwZgGRqCU0TSQ0XnuP8LmAB0JjjPXfrNtpHgti6RtLP5+80079ic/MPzOffDc+kyVL9BRSS9xGwqd/e/hee3L3P3Pd29R/gY6O5K3JJ2Zk+aze173c7aL9YCKGmLSFqK5z7u282sH7Av0CRi+kOJDEyktu197N6s/WItud3U65mIpK9KE7eZXQMMJ0jcLxIM8/kOoMQtdVLkMJw5bXPoOLAjY18bS8suLRl588hUhyciUiPxdHl6MjACWOnuZwEDAR2ySJ1UOgxn4fJCcNiyZgvL3lzG3Afmpjo0EZFaEU/i3uruJUCRmbUEVgHdEhuWSPVEG4bTS5y3rn8rRRGJiNSueIb1nGVmrYB/EFxdvplg8BCROkfDcIpIfRfPxWkXhE/vNrOXgZbuPj+xYYlUTdH2IjIbZ5Kblxs0k5ejYThFpL6I2VRuZkPKP4A2QKPwuUidsHX9Vu494F4+uOUDRtwwYrcxtDUMp4jUJxUdcf+lgnkOHFnLsYhUS5NWTeg8tDNt92lLr2N7AWgYThGpt2IO61mXaFhPKa9oWxHTr5vOQZccRLMOzVIdjohIrarusJ6lK58Rbbo6YJFUWrd4HTP/NpO2e7dl8NmDUx2OiEjSxHNV+dCI500I7umegzpgkRRYvXA17fdtT4e+Hbjoi4to2aVlqkMSEUmqeK4qvyjydXhr2L8SFpFIDAseXcBTpz/FWTPOIu+QPCVtEWmQ4umApbwfgB61HYhILKXXYfQe3Zuf3PQTuuyvwUFEpOGK5xz3VIKryAEygT7A44kMSqTUp898yux7ZnPqc6eSlZPFsN8PS3VIIiIpFc857okRz4uA5e6+orKVzOx+4Dhglbv3C6e1AR4D8oFlwC/dfX0VY5YGpHhHMVvXbWXb+m26elxEhDiayt39LeAzgoFF2hAk73hMBo4qN+0KYJq77w1MC1+L7GL90vUsfm0xAH1/2Zez3ztbSVtEJBRPU/m5wNXAG4ABt5vZ9e5+f0XrufsMM8svN/kEgiFCAR4EpgN/qFLEUu9EDsOZm5dLdststhdu56IvLiKzcSYZmdW5FENEpH6Kp6n8cmCwu68FMLO2wHtAhYk7ho7u/l34fCXQMdaCZjYeGA+Ql5dXjU1JOigdhrN0RK/C5YU0atqIETeMILNxZoqjExGpe+I5lFkLbIp4vSmcViMeXCocs9s2d5/k7gXuXtC+ffuabk7qqGjDcBZtLeKDv32QoohEROq2eI64vwRmmtmzBIn2BGC+mV0K4O5/rcL2vjezPdz9OzPbg2Bsb2nANAyniEjVxHPEvRh4hh+Pjp8FlgItwkdVPAecGT4/MyxLGqDiHcW89vvXaN6pedT5GoZTRCS6eHpOuw7AzJqHrzfHU7CZTSG4EK2dma0ArgFuAh43s3OA5cAvqxe2pLttG7Yx/5/zyTs8jy+mfrFLc7mG4RQRiS2eq8r7Af8kuBUMM1sDnOHun1S0nrufFmOWvpEbsBUzV9Bl/y4069CM8z8+n5y2ObtdVa5hOEVEYovnHPck4FJ3fxPAzIYD/wAOTmBcUg8tfWMpD414iBMfPpEBYwaQ0zYHgP5j+itRi4jEKZ5z3M1KkzaAu08H1BuGxM1Lgssj8ofnc8ydx7DvyfumOCIRkfQVT+JeYmb/z8zyw8f/AEsSHZjUD1++8iWT9pvE1nVbsQxj6PlDaZQdT0OPiIhEE0/iPhtoDzwFPAm0C6eJVCqnbQ6NmjZix+YdqQ5FRKReiHnoY2ZNgF8DewELgN+5+85Yy4uU2rB8A8tnLGfg2IF0LujM2e+ejZmlOiwRkXqhoiPuB4ECgqR9NHBzUiKStPf2jW/zyoRX2LZhG4CStohILaroZOO+7t4fwMzuAz5MTkiSjop3FLN903Zy2uYwcuJIhv1+GE1aNUl1WCIi9U5FibusWdzdi3TUJLG4O48c8wglRSWc+caZZLfIJrtFdqrDEhGplypK3APNbGP43ICm4WsjGCOkZcKjkzonVmcpg88ZTGZWJpahH3giIokUM3G7u8ZUlF1EG4Lz2XOC7ubVgYqISHLEczuYCBB9CM7i7cVMu2paiiISEWl4lLglbhqCU0Qk9ZS4JW4aglNEJPWUuKVSO34Iej0befNIMrN3vfRBQ3CKiCSXErdU6PPnP+dv+X9jzWdr6D+mPyfcdwK53XPBILd7LqMmjdKFaSIiSaTRHqRCnQZ3oseRPWiSG3SmoiE4RURSS0fcspsVM1fw2h9ew91p2aUlJz92cszz2yIiklxK3LKbxa8u5pPHPmHLmi2pDkVERMoxd091DJUqKCjwWbNmpTqMem3zys1sWbOFDv06UFJcwo7NO8qax0VEJLnMbLa7F0SbpyNuwd2ZcvwUnjztSbzEycjMUNIWEamjdHFaA1a8s5iMzAwswzjmjmPIyslSX+MiInVcSo64zewSM/vEzD42sylmpsO7JNtWuI0HDn2A9295H4Au+3ehQ78OKY5KREQqk/TEbWZdgIuBAnfvB2QCpyY7joYuu2U27fu0p1X3VqkORUREqiBVTeWNCIYJ3QnkAN+mKI56L3IYzpZdW9KhfwdGPzCaZh2accIDJ6Q6PBERqaKkH3G7+zfAROAr4Dug0N1fTXYcDUHpMJyFywvBYePXG/nyxS958+o3Ux2aiIhUUyqaylsDJwA9gM5AMzM7Pcpy481slpnNWr16dbLDrBeiDcMJ8OXLX6YgGhERqQ2puDjtJ8BSd1/t7juBp4CDyy/k7pPcvcDdC9q3b5/0IOsDDcMpIlL/pCJxfwUcaGY5ZmbACGBRCuKo13Zu3RlzuE0Nwykikr5ScY57JvAEMAdYEMYwKdlx1GevXv4qDx35EEf87xFk5WTtMk/DcIqIpLeUXFXu7tcA16Ri2w1Bl6FdyMjMoN8p/cjIyCi7qjw3L5cRN4zQ6F4iImlMfZXXAyVFJcz44wza79uevr/sm+pwRESkhtRXeQOw+JXFfPXOV6kOQ0REEkx9lacpd2fBIwvofWJvGjdrzNjXx9K4WeNUhyUiIgmmI+409f2873n6jKeZc+8cACVtEZEGQkfcaab0IrNOgzox7q1x5A3LS3VIIiKSRDriTiOzJ83m9l63s3pR0JNc90O7axhOEZEGRkfcacDdMTN6j+7Nxm820nrP1qkOSUREUkSJu46JHM0rNy+XLgd0gRI4+fGTadahGUdcd0SqQxQRkRRS4q5DSkfzKh0YpHB5IZu+3UTnoZ0p3lFMo2ztLhGRhk7nuOuQaKN5lewsYdM3m5S0RUQEUOKuUzSal4iIVEaJuw4o3lnM7H/MJrebRvMSEZGKKXHXAUteW8Lz45+n9897azQvERGpkBJ3iuzYvIMVM1cAsNfRe3H2u2dz1C1HMWrSKHK754JBbvdcRk0apdG8RESkjEYHS5Enf/Uki19dzITlE9RdqYiI7KKi0cF0qXIS/bD6Bxo1aUR2i2wOv+Zwhl4wVElbRESqRE3lSbKtcBt39r2TN656A4B2+7Qj7xD1My4iIlWjI+4E275pO9ktsmmS24TDrzmc/OH5qQ5JRETSmI64E+iLF7/g1rxbWb0wGBRk/9/sT4e+HVIclYiIpDMl7gTwkuCCv85DO9NrVC+yc7NTHJGIiNQXaiqvofKDguyx3x5QAr986pc0a9+MEx86MdUhiohIPaLEXQOxBgXpNqwbJTtLyGycmeIIRUSkvklJ4jazVsC9QD/AgbPd/f1Eb7f80fGIG0ZUq3OTVR+vIqtZVsxBQTYs3aCkLSIiCZGqI+6/AS+7+8lm1hjISfQGox0dTx0/FYCOgzoClF049u6f36X5Hs0ZOHYgAPcdfB/dDu7GyIkjAZh8+GT6ntpXg4KIiEjSJT1xm1kucBgwDsDddwA7Er3daEfHO7fsZNpV02jauiktu7bktKmnAfDxvz6mY/+OZYm760Fdade7Xdl6Jz58Irl5uXzxwhcULt89SWtQEBERSZSkd3lqZoOAScBCYCAwG/itu/9QbrnxwHiAvLy8/ZYvX16j7V6XcV3QKL9bQHDWjLNo3LwxnQZ1AsDdMbNKyyx/FA/BoCDqX1xERGqioi5PU3E7WCNgCHCXuw8GfgCuKL+Qu09y9wJ3L2jfvn2NNxrrKDg3L5e8Q/LKkjYQV9IG6D+mvwYFERGRpErFOe4VwAp3nxm+foIoibu2jbhhRNSj45oOmdl/TH8lahERSZqkH3G7+0rgazPbJ5w0gqDZPKF0dCwiIvVBqq4qvwh4JLyifAlwVjI2qqNjERFJdylJ3O4+F4h60l1ERERiU1/lIiIiaUSJW0REJI0ocYuIiKQRJW4REZE0kvSe06rDzFYDNes6bVftgDW1WF4qqS51T32pB6gudVV9qUt9qQfUfl26u3vU3sfSInHXNjObFasruXSjutQ99aUeoLrUVfWlLvWlHpDcuqipXEREJI0ocYuIiKSRhpq4J6U6gFqkutQ99aUeoLrUVfWlLvWlHpDEujTIc9wiIiLpqqEecYuIiKQlJW4REZE0kpaJ28y6mdmbZrbQzD4xs9+G09uY2Wtm9kX4t3U43czsNjP70szmm9mQiLLODJf/wszOjLG9qOXWlXqY2SAzez8sY76ZnRJje+PMbLWZzQ0f59ZGPWqzLuG84ogYn4uxvWwzeyxcf6aZ5de1upjZERH1mGtm28xsdJTt1aX90jv8LG03s8vKlXWUmX0W1vOKGNtLyH6prXrEKifK9oabWWHEPrm6NupRm3UJ5y0zswVhjLNibC/m/1pdqYuZ7VPuf2WjmU2Isr2E7Jdq1GNM+F4uMLP3zGxgRFmJ/z9x97R7AHsAQ8LnLYDPgX2BPwNXhNOvAP4UPj8GeAkw4EBgZji9DcGwom2A1uHz1lG2F7XcOlSPXsDe4fPOwHdAqyjbGwfcUZf3SThvcxzbuwC4O3x+KvBYXaxLRJltgHVATh3fLx2AocANwGUR5WQCi4E9gcbAPGDfZO2XWqxH1HKibG848Hxd3ifhvGVAu0q2V+nnsy7UpdxnbSVBByRJ2S/VqMfBhLkCOJofv4uT8n9S6x/KVDyAZ4GfAp8Be0TsiM/C5/cAp0Us/1k4/zTgnojpuyxXfvny5daVekQpZx5hIi83fRwJShC1WRfiS9yvAAeFzxsR9Fhkda0uEdPGA4/EKL/O7JeI5a5l14R3EPBKxOsrgStTtV+qW49Y5USZPpwEJe7arAvxJe64vjdSXZeIeSOBd2PMS8p+ibce4fTWwDfh86T8n6RlU3mksIlhMDAT6Oju34WzVgIdw+ddgK8jVlsRTos1vbxY5daaGtYjspz9CX7pLY6xqZPCJp4nzKxb7US/q1qoSxMzm2VmH1iUpuXy67t7EVAItK2tOpSqrf1C8Kt6SgWbqiv7JZZ4/1cSvl9qWI9Y5URzkJnNM7OXzKxvdeOtQgzVqYsDr5rZbDMbH2OZePddjdTWfqHy/5WE7pdq1OMcghYNSNL/SVonbjNrDjwJTHD3jZHzPPgpU+v3uiWi3Nqqh5ntAfwTOMvdS6IsMhXId/cBwGvAgzUKPHoMtVGX7h50Hfgr4FYz61nbccajlvdLf4Jf2dGky35JuVrcJzHLCc0h+BwOBG4HnqlR4FWMoQp1OcTdhxA01/7GzA6r7TjjUYv7pTFwPPDvGIskdL9UtR5mdgRB4v5DbcZRmbRN3GaWRfAGP+LuT4WTvw+/JEu/LFeF078BIo9iuobTYk0vL1a5daUemFlL4AXgKnf/INq23H2tu28PX94L7Fdb9ajNurh76d8lwHSCX7/lla1vZo2AXGBtXatL6JfA0+6+M9q26th+iSXe/5WE7Zdaqkescnbh7hvdfXP4/EUgy8za1UI1KoqhynWJ+F9ZBTwN7B9lsXj3XbXUVl1CRwNz3P37aDMTuV+qWg8zG0Dw/3qCu5d+xpPyf5KWidvMDLgPWOTuf42Y9RxwZvj8TILzFKXTz7DAgUBh2PzxCjDSzFqHVwuOJPpRUaxy60Q9wl+pTwMPufsTFWxvj4iXxwOLaqMeYdm1VZfWZpYdltkOGAYsjLLJyHJPBt4IfxHXmbpErHcaFTT91bH9EstHwN5m1iP8vJ0allFeQvZLbdWjgnLKL9cpXLb09FMGtfcDpLbq0szMWpQ+J/j++jjKopV9PqutFj9fpSr7X0nIfqlqPcwsD3gKGOvun0csn5z/k3hPhtelB3AIQZPFfGBu+DiG4BzBNOAL4HWgTbi8AX8nOO+7ACiIKOts4MvwcVbE9HtLl4tVbl2pB3A6sDOijLnAoHDe9cDx4fP/Az4huHjtTaB3XdsnBFdrLghjXACcE7GNyLo0IWhO+xL4ENizrtUlnJdP8Os6o9w26up+6URwXm4jsCF83jKcdwzB1baLCVp2krZfaqsescoJ1/k18Ovw+YUR++QD4OC6tk8IrlyeFz4+KbdPIusS8/NZV+oSzmtGkIRzy20j4fulGvW4F1gfseysiLIS/n+iLk9FRETSSFo2lYuIiDRUStwiIiJpRIlbREQkjShxi4iIpBElbhERkTSixC2SxsyslZldEPG6s5nFvJe/htsabbUwGpOZ9TezybUQkkiDpNvBRNKYBf0qP+/u/ZKwrfcI7kNdE+fyjTzohznavNeBs939q9qMUaQh0BG3SHq7CehpwdjEN5tZvpl9DGXjfD9jwTjCy8zsQjO71Mz+Y8EALm3C5Xqa2csWDFTxtpn1Lr8RM+sFbHf3NWbWwsyWhl1EYmYtS1+b2XQzu9WCsaF/a2a/MLOPLRgUYkZEkVMJepUSkSpS4hZJb1cAi919kLtfHmV+P+Dn/DgG8hZ3Hwy8D5wRLjMJuMjd9wMuA+6MUs4wggEecPdNBH3IHxvOOxV4yn/si72xuxe4+1+Aq4GfeTAoxPER5c0CDq1GfUUaPCVukfrtTXff5O6rCYYOnBpOXwDkWzAa0sHAv81sLsHYzXtEKWcPYHXE63uBs8LnZwEPRMx7LOL5u8BkMzsPyIyYvgroXL0qiTRsjVIdgIgk1PaI5yURr0sI/v8zgA3uPqiScrYSjGAEgLu/GzbLDwcy3T1ycIsfIpb7tZkdQHB0PtvM9vNgJKUmYZkiUkU64hZJb5uAFtVd2YMxh5ea2S8gGCXJzAZGWXQRsFe5aQ8Bj7Lr0fYuzKynu89096sJjthLhzzsRfSRrESkEkrcImksPHp9N7wA7OZqFjMGOMfMSkeZOiHKMjOAwaVDKoYeAVpTwTCMwM1mtiC8YO49glGdAI4gGD9eRKpIt4OJSFzM7G/AVHd/PXx9MnCCu4+tYjnZwFvAIbFuFxOR2HSOW0TidSNwAICZ3Q4cTTD2cFXlAVcoaYtUj464RURE0ojOcYuIiKQRJW4REZE0osQtIiKSRpS4RURE0ogSt4iISBr5/4U0S73v4WUwAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "fig = plt.figure(figsize=(8,4))\n", "plt.plot(t,w,'o:',color='purple',label='Adams-Bashforth')\n", "plt.title('Population Equation with seasonal oscilation')\n", "plt.xlabel('time (yrs)')\n", "plt.ylabel('Population in Billions')\n", "plt.legend(loc='best')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": { "id": "6AXoYl16aNb9" }, "source": [ "## Table\n", "The table below shows the time and the numerical approximation, $w$, for the non-linear population equation with oscilations:" ] }, { "cell_type": "code", "execution_count": 28, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 696 }, "id": "1WTdnm0PaNb-", "outputId": "7887162a-a28a-4812-cf80-51dd5f83e8fd" }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
time t_iAdams Approx w
02000.06.000000
12001.06.110000
22002.06.963018
32003.07.900330
42004.08.880317
52005.09.883555
62006.010.889620
72007.011.877816
82008.012.828881
92009.013.726473
102010.014.558187
112011.015.315964
122012.015.995957
132013.016.597982
142014.017.124739
152015.017.580977
162016.017.972719
172017.018.306611
182018.018.589435
192019.018.827758
202020.019.027711
\n", "
" ], "text/plain": [ " time t_i Adams Approx w\n", "0 2000.0 6.000000\n", "1 2001.0 6.110000\n", "2 2002.0 6.963018\n", "3 2003.0 7.900330\n", "4 2004.0 8.880317\n", "5 2005.0 9.883555\n", "6 2006.0 10.889620\n", "7 2007.0 11.877816\n", "8 2008.0 12.828881\n", "9 2009.0 13.726473\n", "10 2010.0 14.558187\n", "11 2011.0 15.315964\n", "12 2012.0 15.995957\n", "13 2013.0 16.597982\n", "14 2014.0 17.124739\n", "15 2015.0 17.580977\n", "16 2016.0 17.972719\n", "17 2017.0 18.306611\n", "18 2018.0 18.589435\n", "19 2019.0 18.827758\n", "20 2020.0 19.027711" ] }, "execution_count": 28, "metadata": {}, "output_type": "execute_result" } ], "source": [ "d = {'time t_i': t, 'Adams Approx w': w}\n", "df = pd.DataFrame(data=d)\n", "df" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "id": "HgiQsWs2aNb_" }, "outputs": [], "source": [] } ], "metadata": { "colab": { "include_colab_link": true, "name": "402_Adams Bashforth Population Equations.ipynb", "provenance": [] }, "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.11.4" } }, "nbformat": 4, "nbformat_minor": 1 }